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ABSTRACT 


Active  damping  of  modal  oscillation  is  critical  to  the 
success  of  future  versions  of  Space  Station  Freedom. 
Vibratory  motion  may  be  induced  by  external  disturbances  such 
as  solar  and  gravity  gradient  torques,  extra  vehicular  and 
experimental  activity,  aerodynamic  forces,  the  earth's 
magnetic  field,  and  space  shuttle  docking.  Linear  proof  mass 
actuators  can  provide  control  on  the  space  station  to  achieve 
this  damping  effect.  Two  control  algorithms,  Linear  Quadratic 
Gaussian  control  and  H®  control  are  applied  to  a  model  of 
Space  Station  Freedom.  The  results  compare  the  robustness, 
stability,  and  performance  of  the  Space  Station  under  the 
effects  of  each  of  the  two  control  algorithms. 
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I .  INTRODUCTION 

A.  OVERVIEW 

A  long  time  quest  for  scientists  and  engineers  has  been 
the  development  of  a  permanent  space  station.  This  pursuit 
has  resulted  in  the  ongoing  development  of  Space  Station 
Freedom  which  is  intended  to  be  a  permanently-manned  orbiting 
base  by  the  end  of  the  century.  The  identification  of  the 
placement  of  actuators  and  sensors,  diminishing  effects  of 
extra  vehicular  activity  and  zero  gravity  on  crew  members, 
developing  of  solar  panels  for  energy,  and  damping  structural 
vibrations  are  some  of  the  tasks  and  obstacles  involved  in  the 
design  of  Space  Station  Freedom.  Of  particular  interest  is 
the  control  of  vibratory  motion  on  the  space  station. 
Classical  control  techniques  have  been  employed  in  the  past  to 
solve  the  vibration  problem.  This  approach  alone  is  limited 
due  to  its  inability  to  address  the  issues  of  robustness, 
sensitivity,  and  fault  tolerance  in  a  multivariable  setting. 
Modern  advances  in  control  systems  and  computer  technology 
have  equipped  designers  with  the  necessary  tools  to  overcome 
these  limitations. 

The  present  study  is  aimed  at  developing  a  control 
algorithm  for  damping  structural  vibrations  in  the  presence  of 
parameter  variations  and  unmodelled  dynamics.  In  order  to 
achieve  this  goal,  H°o  control  and  Linear  Quadratic  Gaussian 
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control  will  be  evaluated  and  applied  to  a  model  of  Space 
Station  Freedom. 

B.  SPACE  STATION  FREEDOM  AND  CONTROL  SYSTEMS 

Space  travel  has  been  envisioned  by  astronomers  and 
scientist  since  approximately  300  A.D.  [Ref.  1].  Early 
references  to  space  travel  were  only  stories  with  no 
scientific  foundation.  Later  in  the  1700's,  Galileo  first 
viewed  the  stars  and  moons  when  he  envisioned  these  heavenly 
bodies  as  places  for  man  to  visit.  About  1860,  Jules  Verne 
detailed  the  venture  of  three  crew  members  who  soared  in  a 
spacecraft  to  the  moon  in  De  La  Terre  a'  la  Lune  (From  the 
Earth  to  the  Moon) .  The  ideas  of  space  travel  persevered  into 
the  1900's.  After  World  War  II,  the  space  station  concept  and 
serious  proposals  for  a  manned  spacecraft  originated.  Manned 
spaceflight  became  a  national  goal  after  President  Kennedy's 
pronouncement  to  send  a  man  to  the  moon.  Hence,  the  Mercury, 
Gemini,  and  Apollo  space  programs  resulted. 

The  manned  lunar  landing  was  a  prominent  event  in  the 
space  program.  The  diversion  of  attention  to  the  lunar 
landing  missions  led  to  the  decline  of  interest  toward 
building  a  permanent  space  station.  Later,  mandated  budgetary 
constraints  imposed  by  Congress  did  not  improve  the  situation. 
In  1982,  Space  Station  Freedom  became  a  national  goal  [Refs. 
2  and  3].  Development  of  Space  Station  Freedom  was  divided 
into  eight  phases.  Control  and  stabilization  were  two  phases 
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that  focused  on  limiting  the  effects  of  structural  vibratory 
motion. 

The  beginnings  of  control  theory  has  been  traced  from 
Huygins  in  1673  to  Maxwell  in  1868  [Ref.  4],  Maxwell  was  the 
first  to  discuss  stability  of  feedback  control  systems  whereas 
Bode,  Routh,  and  Nyquist  later  addressed  this  issue.  Feedback 
control  became  an  important  topic  at  approximately  the  same 
time  the  space  concept  emerged.  Control  design  was  prominent 
due  to  the  construction  of  rockets,  missiles,  and 
communication  systems  developed  during  WWII.  This  effort 
undoubtedly  impacted  the  space  design  process.  With  the  need 
for  space  control  technology,  Kalman,  Pontryagin,  and  Bellman 
investigated  at  optimal  control  and  state  space  modelling 
[Refs.  5  and  6]  while  Joseph,  Tou  and  Simon  developed 
compensators  that  estimated  the  states  of  a  system  in  the 
presence  of  stochastic  noise  [Ref.  7]. 

Most  of  the  systems  developed  during  this  time  were 
single-input/single-output  (SISO) .  However,  the  requirement 
for  controlling  systems  with  multiple  inputs  and  multiple 
outputs  (MIMO)  was  becoming  important.  Zames  [Ref.  8] 
introduced  a  feedback  control  design  which  addressed 
multivariable  systems  affected  by  external  perturbations. 
Classical  and  feedback  control  techniques  in  the  frequency 
domain  were  extended  to  MIMO  systems  with  the  advent  of  the  H°o 
(H-infinity)  design  methods.  The  H«  design  methodology  allows 
the  designer  to  directly  consider  the  contradictory 
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requirements  of  system  performance,  sensitivity  reduction, 
robustness  and  disturbance  attenuation  in  multivariable 
control  systems.  The  H°o  controller,  is  consequently,  a 
natural  choice  for  suppression  of  vibrations  due  to 
disturbances  and  unmodelled  dynamics  aboard  Space  Station 
Freedom. 

C.  VIBRATION  SUPPRESSION  ON  SPACE  STATION  FREEDOM 

Active  damping  of  modal  oscillation  is  critical  to  the 
success  of  future  versions  of  Space  Station  Freedom  [Ref.  9]. 
The  vibratory  motion  may  be  induced  by  external  disturbances 
such  as  solar  and  gravity  gradient  torques,  extra  vehicular 
and  experimental  activity,  aerodynamic  forces,  the  earth's 
magnetic  field,  and  space  shuttle  docking.  Damping  this 
vibratory  motion  reduces  the  disruption  of  onboard  experiments 
and  communication  and  remote  sensor  pointing  errors.  Linear 
proof  mass  actuators  can  provide  control  on  the  space  station 
to  achieve  this  damping  effect  [Ref.  10]. 

Currently,  there  is  little  documented  research  that 
provides  an  analysis  and  comparison  of  different  control 
algorithms.  It  is,  therefore,  the  intent  of  this  research  to 
apply  two  control  algorithms,  Linear  Quadratic  Gaussian 
control  and  H°o  control  to  a  model  of  Space  Station  Freedom. 
The  results  will  compare  the  robustness,  stability,  and 
performance  of  Space  Station  Freedom  under  the  effects  of  each 
of  the  two  control  algorithms. 
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D.  ORGANIZATION  OF  THESIS 


Linear  Quadratic  Gaussian  control  and  H°o  control  are  two 
algorithms  used  to  solve  the  vibration  problem.  Chapter  II 
provides  the  mathematical  model  of  a  space  station.  The  data 
used  in  this  study  are  furnished  by  the  McDonnell  Douglas 
Astronautics  Company.  The  two  control  algorithms  are 
presented  in  Chapter  II.  Because  of  the  complexity  and  fngh 
order  of  the  model,  reduced-order  controllers  are  also 
presented  in  Chapter  II.  Chapter  III  examines  the  application 
of  the  two  control  systems  to  the  model  and  a  comparison  of 
the  controllers.  Conclusions  and  recommendations  for  future 
study  are  presented  in  Chapter  IV. 
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II.  CONTROL  SYSTEM  DESIGNS 


The  space  station  exhibits  low  frequency,  lightly  damped 
vibrations  due  to  its  size,  construction  and  composition. 
Figure  1  is  a  representation  of  a  dual-keel  space  station 
provided  courtesy  of  the  McDonnell  Douglas  Astronautics 
Company . 


WORK  PACKAGE-2  ELEMENTS  AND  SYSTEMS*^™ 


Figure  l  a  Mathematical  Model  of  Space  Station  Freedom 
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A  mathematical  model  is  required  to  analyze  the  motion  of 
this  space  station.  The  modal  equation  can  be  found  by 
applying  finite  element  analysis  to  the  space  station 
dynamics.  The  motion  is  then  a  linear  combination  of  tne 
natural  modes.  The  solution  of  the  modal  equation  is  of  the 
form  [Ref .  11 ] : 

00 

git)  =  ( t)  (2-1) 

j  =i 


where  q(t)  is  the  displacement  of  the  structure,  is  the 
modal  vector  or  natural  mode  shape,  and  77 i  ( t )  is  the  modal 
amplitude  of  the  ith  mode  at  time  t.  The  modal  amplitudes  can 
be  found  by  solving  the  uncoupled  second-order  differential 
equations : 

iii  ( t)  +  d  wx  f)i  +  t)  ( t)  =  $i(p)  •  fp(t.)  (2-2) 

where  d  is  a  scalar  structural  damping  coefficient,  is  the 
ith  natural  modal  frequency  and  f  is  an  external  forcing 
function.  Solving  for  the  highest  order  derivative  in  Equation 
2-3  to  obtain  a  state  space  representation  of  the  modal  matrix 
gives : 


=  -dw\  Tji(t)  + 

=  Qtip)  •  fp(  t)  i  =  1,2,  .  .  .10 


(2-3) 
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and  let 


( t)  =  ( t) 

jq  ( t)  =  Xj  =  i| j  ( t) 

x^(t)  =  -dwix2-w?ix1  +  QI(p)  ■  fx{t)  =  i\(t) 


(2-4) 


where 


fp{t)  =  u(  t) 

^(p)  =  b^p)  . 


(2-5) 


The  matrix  state  equation  is: 


0  1 

0 

.’ii. 

-w'\  -dwi 

V 

bAp) 

A  measurement,  taken  at  a  point  p,  is  given  by: 


y(t)  =  C  x(t) 


’ll 


=  [<I>1  (p)  0  .  .  .  3>10  (p)  0] 


Tlio 


(2-7) 


The  data  provided  100  modes  of  vibration  and  their 
corresponding  natural  frequencies  [Ref.  12].  In  this  study, 
however,  only  the  first  ten  modes  are  considered.  Tables  1 
and  2  present  these  ten  modes  and  their  natural  frequencies, 
respectively.  Vibration  may  occur  in  the  translational  or 
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rotational  direction.  Furthermore,  node  23,  the  shuttle 
docking  point  is  analyzed. 


TABLE  1 

LIST  OF  TRANSLATIONAL  MODAL  VECTORS 


MODE 

NODE  23 

X 

Y 

Z 

1 

2 . 7006e-02 

-2 . 7171e-13 

7 . 8058e-12 

2 

-3 . 2282e-13 

2 . 7006e-02 

4 . 700Ge-13 

3 

1 . 7084e-12 

1 . 5979e-13 

2 . 7006e-02 

4 

2 . 6363e-13 

3 . 2302e-03 

3 . 3864e-03 

5 

-4 . 9299e-03 

2 . 1923e-05 

2 . 52  67e-02 

6 

-4 . 2493e-03 

2 . 0629e-02 

1. 3326e-04 

7 

-1. 4465e-03 

-5 . 9130e-04 

-7 . 28 15e-03 

8 

6 . 4  3  08e-03 

3 . 7219e-04 

-1 . 9019e-03 

9 

-7 . 0786e-05 

3 . 6509e-04 

-1 . 2529e-04 

10 

-1 . 07  88e-03 

2 . 6420e-04 

2 . 64  2  0e-04 
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TABLE  2 

NATURAL  FREQUENCIES 


Modes 

1 

0.32374935 

2 

0.34734376 

3 

0.37886664 

4 

0.38108557 

5 

0.39618959 

6 

0.40531164 

7 

0.40761414 

8 

0.41520832 

9 

0.43176959 

10 

0.44403511 

A.  CONTROL  ALGORITHMS 

There  are  a  number  of  linear  control  algorithms  available 
to  a  designer.  The  Linear  Quadratic  (LQ)  regulator,  Kalman 
filter,  Linear  Quadratic  Gaussian  (LQG)  control,  and  H=» 
control  are  applied  and  analyzed  in  this  thesis.  The  Loop 
Transfer  Recovery  ( LTR)  technique  is  also  presented. 

1.  The  Linear  Quadratic  Regulator 

The  LQ  regulator  is  an  algorithm  that  results  in  an 
optimum  controller  for  a  given  system.  A  performance  measure 
is  minimized  given  that  the  states  of  the  system  are  fully 
accessible.  Kirk  [Ref.  13  ]  states  that  the  objective  is  to 
determine  the  control  signals  that  will  cause  a  process  to 
satisfy  the  physical  constraints  and  at  the  same  time  minimize 
(or  maximize)  some  performance  measure.  The  selection  of  the 
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performance  measure  is  based  on  the  objectives  and  imagination 
of  the  designer. 

In  addition  to  the  performance  measure,  the 
formulation  of  the  LQ  regulator  requires  a  mathematical  model 
of  the  space  station.  The  mathematical  model  utilizes  the 
linear,  time-invariant  state  equation: 

x{k+ 1)  =  $  x(k)  +  r  u(k )  (2-8) 

where  *  is  an  n  by  n  matrix,  r  an  n  by  r  matrix,  x(k)  an  n- 
dimensional  state  vector.  u(k)  a  r-dimensional  control  vector 
which  minimizes  the  performance  measure: 

N-l  ' 

J  =  £  f xT(k)  Q  x(k)  +  uT(k)  R  u(k )]  (2-9) 

*r  =  C 


in  which  Q  is  the  square,  symmetric  state  weighting  matrix  and 
R  is  the  control  weighting  matrix.  To  guarantee  a  unique 
solution,  Q  must  be  positive  semi-definite  and  R  positive 
definite.  The  Q  matrix  provides  a  penalty  for  deviation  from 
equilibrium  while  R  provides  a  penalty  for  using  control. 
Consequently,  increasing  Q  increases  the  penalty  on  the  state 
vector,  while  increasing  R  increase  the  cost  of  control 
applied  to  the  system.  The  performance  function  is  minimized 
when  the  linear  optimal  control  law,  u(k)  =  -Kx(k)  is 
implemented.  The  optimal  gain  matrix  K  is  the  solution  of  the 
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algebraic  Ricatti  equation  that  drives  the  states  to  zero. 
Figure  2  is  a  block  diagram  of  the  system: 


U(k) 

A 

i 


i 


PLANT 


Y(k)  . 

~7— > 

i 

I 

i 

i 


Figure  2  The  Linear  Quadratic  Regulator 

The  LQ  regulator  has  good  stability  and  robustness 
properties,  but  assumes  that  no  noise  exists  in  the  system  and 
all  states  are  available  for  measurement  and  feedback.  This 
assumption  poses  a  major  drawback  for  the  LQ  regulator  and  may 
be  overcome  by  applying  the  Kalman  filter. 

2.  The  Kalman  Filter 

The  Kalman  Filter  is  an  observer  that  provides  optimal 
state  estimation  in  the  presence  of  noise.  It  is  a  recursive 
algorithm  that  provides  the  "best"  estimate  of  the  states  of 
a  linear  system  given  only  the  input  and  output  of  the  system. 
The  linear,  time-invariant  system  whose  state  is  to  be 
estimated  is: 
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(2-10) 


jr(k+l)  =  O  x(Jc)  +  I\  u(k)  +  r2  w{k) 
y(k)  =  C  x(k)  +  v(k) 

where  #  is  an  n  by  n  matrix,  and  r2  are  n  by  r  matrices, 
u(k)  is  a  known  input,  x(k)  is  a  state  matrix,  w(k)  is  random 
plant  driving  noise  and  v(k)  is  random  measurement  noise.  Both 
w(k)  and  v(k)  are  white  noise  vectors  with  zero  mean,  i.e., 
E[w]  -  E[v]  =  0.  The  covariance  matrix  for  the  plant  driving 
noise  is  E[wwT]  =  W  and  the  measurement  noise  covariance 
matrix  is  E[wT]  =  V.  The  plant  driving  noise  and  measurement 
noise  are  independent  and  uncorrelated. 

The  optimal  state  estimate: 

£(k+l\k*l)  =x(k*l\k)  +  G[y(k+1)  -  y{k+l  \  k)  ]  (2-11) 


where  Jt(k+l]k)  denotes  the  estimate  of  the  state  at  time  k+1 
given  measurements  through  time  k  and  y(k+l|k)  is  the  estimate 
of  the  measurement  at  time  k+1  given  measurement  through  time 
k.  A  summary  of  the  other  equations  that  represent  a  steady- 
state  Kalman  filter  are: 


jt(k  +  l\k)  =  $x(k\k)  +  I\  u(k)  +  Ta  w(k) 
?(k+ 1 ! k)  =  Cx(k+1 i k) 


(2-12) 


A  general  block  diagram  of  a  Kalman  filter  is  given  in  Figure 
3  [Ref.  14]: 
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Figure  3  The  Kalman  Filter 

The  Kalman  filter  is  an  optimum  process  provided  the 
stochastic  noise  processes  are  white  and  Gaussian  [Ref.  12]. 
3.  Linear  Quadratic  Gaussian  Control 

The  Linear  Quadratic  Gaussian  (LQG)  control  is  a 
combination  of  the  LQ  regulator  and  the  Kalman  filter  designed 
in  separate  stages.  The  results  derived  separately  for  the 
optimal  control  and  estimation  problem  are  still  valid.  A 
typical  block,  diagram  of  the  system  is  depicted  in  Figure  4. 
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Figure  4  The  Linear  Quadratic  Gaussian  Controller 

The  LQG  control  system  utilizes  the  full-state  feedback  of  the 
LQ  regulator  applied  to  the  estimated  states  from  the  Kalman 
filter.  The  robustness  properties  of  the  feedback  system 
diminishes  when  the  Kalman  filter  is  implemented  to  estimate 
the  states  [Ref.  15].  Two  criteria  that  can  be  used  to  assess 
the  robustness  of  the  system  are  the  phase  and  gain  margins. 
The  gain  margin  is  the  amount  of  additional  gain  that  can  be 
added  to  the  system  without  causing  instability.  The  phase 
margin  is  the  amount  of  additional  phase  delay  in  the  plant 
and/or  compensator  that  can  be  added  before  generating 
instability.  In  order  to  improve  the  robustness  as  defined  by 
the  phase  and  gain  margins,  loop  transfer  recovery  techniques 
may  be  applied  [Ref.  16]. 

Loop  transfer  recovery  ( LTR)  is  a  technigue  to  recover 
the  robustness  properties  lost  due  to  using  state  estimates 
acquired  by  the  Kalman  filter.  In  order  to  apply  this 
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algorithm,  the  system  must  be  causal  and  have  minimum  phase, 
i.e.,  all  the  zeros  strictly  lie  inside  the  unit  circle. 
Fictitious  noise,  describing  the  unmodelled  dynamics,  is 
implanted  at  the  plant  input  thereby  adjusting  the  estimator 
design.  Consequently,  the  LQG  controller  becomes  more  robust 
to  gain  and  phase  changes  at  the  plant  input  as  the  noise 
increases  [Ref.  17].  However,  robustness  is  inversely  related 
to  the  performance.  The  system  becomes  suboptimal  with  the 
addition  of  the  fictitious  noise.  This  forces  the  designer  to 
form  a  compromise  between  robustness  or  performance. 

4 .  Hoo  Control 

Since  no  design  model  can  approximate  a  physical  plant 
perfectly,  modelling  errors  will  always  exist.  The  LQ 
regulator,  LQG  control,  and  the  Kalman  filter  possess 
reasonable  tolerance  to  modelling  errors.  However,  these 
control  techniques  lack  the  provisions  to  directly  design  for 
robustness  to  unmodelled  dynamics.  The  Hoo  controller  is  an 
optimal  method  that  is  capable  of  addressing  these  factors 
using  frequency  domain  techniques.  The  Hoo  performance 
function  puts  limits  on  the  maximum  value  of  the  disturbance 
frequency  response. 

The  Hoo  controller  is  similar  in  structure  to  that 
found  in  the  classical  control  design  as  illustrated  in  Figure 
5.  G(s)  is  the  open  loop  transfer  function,  K(s)  is  the 
controller,  U  is  the  known  input,  and  Y  is  the  error. 
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>  G(s) 
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Figure  5  Closed  Loop  System  with  H°o  Controller 

The  objective  of  H<»  is  to  achieve  the  desired  robustness  and 
performance  response  in  the  presence  of  unmodelled  dynamics. 
An  analysis  of  the  response  of  the  multivariable  system  in 
Figure  5  may  be  ascertained  by  an  inspection  of  the  closed 
loop  transfer  function  matrix  singular  values  [Ref.  18]  which 
will  be  denoted  as  o[- ].  H«  employs  these  singular  values  to 
generate  the  performance  measure.  The  performance  measure  is 
the  maximum  singular  value  of  the  closed  loop  transfer 
function  over  a  frequency  range,  i.e.,  minimizing  the  H°°  norm. 
The  singular  values  of  the  closed  loop  system  can  be  used  to 
quantify  its  stability  margins  [Ref  18].  The  conceptual  H°° 
design  procedure  consists  of  four  basic  operations  which  are 
described  below. 

First,  the  design  problem  is  formulated  as  a  "standard 
problem"  as  depicted  in  Figure  5  where  the  objective  is  to 
design  a  controller  such  that  the  closed  loop  system  is 


17 


internally  stable.  Minimizing  the  H°o  norm  of  the  transfer 
function  from  the  disturbance  to  the  output  and  the 
measurement  noise  to  the  output,  is  also  an  objective. 
Secondly,  an  analysis  to  determine  the  sensitivity  of  the 
system  is  required.  The  sensitivity,  S(s),  is  a  measure  of 
how  much  the  closed  loop  transfer  function  changes  with  a 
small  change  in  G(s)  and  is  defined  as  [I  +  G(s)K(s)]  1  = 
¥(5)711,(5).  The  complementary  sensitivity,  T(s)  = 

[I  +  G(s)K(s)]  1  G(s)K(s)  =  -  Y(s)/U2(s),  where  U,(s)  and  U2(s) 
are  the  disturbance  and  measurement  noises,  respectively.  A 
graphical  representation  is  illustrated  in  Figure  6. 


g(s) 


Ul(s) 

V  Y(S) 

>  - > 


K(s)  <- 


V  U2(s) 

<£< - 


Figure  6  Closed  Loop  System  with  Disturbance  and 
Measurement  Noise 

The  sensitivity  determines  the  disturbance  attenuation.  This 
attenuation  performance  specification  may  be  written  in  the 
real  frequency  domain  as  a(S(ju))  <  j  W,  ’ 1  ( j  o )  |  where  |w1'1(jw)| 
is  the  decired  disturbance  attenuation  or  performance  factor 
and  a  denotes  the  largest  singular  value.  Errors  from 
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multiplicative  perturbations  may  affect  the  robustness  whereas 
the  specif  ication  a  (T(  jo)  )  <  |  W3  1  ( j  )  [  where  |w31(jo)|  limits 
the  effects  of  measurement  noise  and  provides  robustness  in 
the  presence  of  these  perturbations.  Figure  7  depicts  the 
closed  loop  system  with  added  weighting  functions,  W^s)  and 
W3  ( s )  . 


Uc  (s) 


>  G(s) 


™§*W3(s) 


Z1  (s) 


v  u(s> 

;  < - 


K  (s) 


Y2(s) 


->  Wl(s) 


Z2(s) 


Figure  7  Augmented  System  with  Weighting  Functions 


The  input/ output  behavior  of  the  system  is  given  by: 


y  ( s) 


-G(s)  K  is)  [I  +  G(s)iC(s)]':l/r?(s)  +  [I  +  G(s)  K(s)  ]  _1  d(s) 
-T(s)  m(s)  +  S{s)  d{s) 

(2-13) 


Third,  an  augmented  system  containing  the  plant  and 
the  weighting  functions  are  produced.  This  augmented  system 
can  sometimes  be  of  high-dimension  and  reduced  order  modelling 
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would  be  warranted  (to  be  discussed  in  section  II. C). 

Fourth,  the  controller  is  designed  for  the  system. 
A  controller  is  selected  that  stabilizes  the  system  and  has  an 
H°°  norm  less  than  one  for  the  closed  loop  system.  This 
restriction  on  the  H®  norm  guarantees  that  the  sensitivity  and 
complementary  sensitivity  meet  the  given  specifications.  The 
resulting  closed  loop  system  is  illustrated  in  Figure  8. 


U(s) 


[Zl(sX  Z2 


Uc(s) 


P(S) 
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- 

J 
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Figure  8  The  H  ®  Controller 


B.  REDUCED  ORDER  MODELLING 

Most  mathematical  models  of  space  structures  are  of  high 
order  which  requires  a  large  computational  effort,  CPU  time, 
and  complex  hardware.  An  alternative  is  to  reduce  the  order 
of  the  model,  thereby  alleviating  some  of  the  complexities  of 
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the  design  problem.  The  control  design  is  then  synthesized 
using  the  reduced-order  model.  The  reduced  order  design 
methods,  however,  do  not  always  produce  feedback  designs  that 
remain  stable  in  the  presence  of  unmodelled  dynamics.  Three 
approaches  to  model  reduction  are  (1)  approximation  of  the 
plant  with  a  low  order  model,  (2)  balance  realizations  by 
organization  of  state  by  order  of  controllability  and  by  (3) 
truncation  of  frequencies.  There  is  no  simple  guideline  to 
determine  which  model  reduction  technique  will  produce  the 
best  result  or  what  degree  of  approximation  is  necessary  to  be 
applied  to  this  model.  The  reduced  order  controller  for  this 
thesis,  however,  will  be  obtained  by  truncation  of 
frequencies . 
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III.  ANALYSIS  AND  SIMULATION 


This  chapter  is  concerned  with  the  design  of  the  LQG  and 
H°o  controllers  for  the  model  discussed  in  Chapter  II.  The 
response  of  the  model  due  to  external  forces  when  both  control 
algorithms  are  applied  is  also  presented.  Two  techniques, 
full  and  reduced  order  LQG  and  H«  controllers,  were  employed 
for  analysis  and  simulations.  Finally,  a  comparison  of  both 
control  algorithms  is  provided. 

A.  FULL  ORDER  CONTROLLERS 
1.  Performance  Measure 

The  objective  for  implementing  the  LQG  and  H<»  control 
is  to  find  a  stable  controller  that  minimizes  the  performance 
measure.  The  total  vibrational  energy  plus  a  control  energy 
term  represent  the  performance  measure.  The  potential  and 
kinetic  energy  equations  written  in  terms  of  the  modal 
amplitudes  and  modal  velocities  are: 

i  10 

P.E.(k)  =  ~Y,w^{k)  (3.1) 

and 

i  10 

K.E.  ( k )  =  -±£  f^dc)  .  (3.2) 

^  2-1 
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The  total  energy  consists  of  the  potential  and  kinetic 
energies  of  the  structure  at  any  point  in  time  and  can  be 
written  in  terms  of  modal  state  vectors  as: 

T.E.  (k)  =  P.E.  (k)  +  K.E.  (k) 

=  -|£  hi  (k)  fk  (k)  ]  Q± 


rii  (k) 
*n  j  (k) 


(3.3) 


where  Q,  is  the  state  weighting  matrix  defined  as: 


Ci 


o>1  0 
0  1 


(3.4) 


Adding  the  control  energy,  the  performance  measure  Equation 
2.10  may  be  written  as: 

10 

J  =  x/cori  +  utRu  (3.5) 

i  =  1 


where  R  is  a  positive  definite  matrix. 

The  performance  measure  of  the  system  may  be 
determined  from  the  system  response  to  a  disturbance.  The 
disturbance  of  interest  is  an  impulse  represented  as  <5  [  •  ]  , 
called  the  delta  or  Dirac  function. 

The  impulse  can  be  used  to  generate  the  initial 
conditions  for  the  system  which  can  be  found  following  the 
derivation  in  Reference  19.  The  state  equation  with  the 
impulse  input  is  given  as: 
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x{t)  =  Ax(t)  +  B5  ( t)  . 


(3.6) 


Integrating  both  sides  of  Eguation  3.6  using  -e  and  e  as  the 
lower  and  upper  limits  respectively: 

x(e)  -  x(-e)  =  B  +  j  * Ax{t)  dt .  (3.7) 

Taking  the  limit  as  e  approaches  0  to  obtain  x(0+),  where  0* 
indicate  the  initial  condition  an  arbitrarily  short  time  after 
the  impulse  occurs,  yields: 


x(0+)  =  B. 


(3.8) 


2.  Open  Loop  System 

The  open  loop  system  is  analyzed  to  assist  in 
understanding  the  closed  loop  system.  The  stability  of  the 
open  loop  system  is  determined  from  the  eigenvalues  of  the 
plant  matrix.  Appendix  A  lists  these  negative,  complex 
conjugate  poles.  The  Bode  plot  of  the  plant  also  illustrates 
the  stability  of  the  open  loop  system  as  depicted  in  Figure  9. 
The  zeros  of  the  plant  are  also  listed  in  Appendix  A.  Note 
that  the  plant  is  a  minimum  phase  system. 
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Figure  9  Open  Loop  Plant  Bode  Plot 


3.  Linear  Quadratic  Gaussian  Controller 

a.  Simulation  of  the  Linear  Quadratic  Regulator 

The  LQG  controller  is  based  on  the  separation 
principle  in  which  the  LQ  regulator  and  the  Kalman  filter  are 
constructed  separately.  The  steady-state  LQG  controller 
utilizes  two  gain  matrices  by  Equation  3.9 


x(k+l) 

.  f* 
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x(k) 
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..*(*+ 1). 

[gc 

F-TK 
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+• 
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out 


( k ) 


[0  K] 


x(k) 

*(k). 


u(k) 


(3.9) 


where  G  and  K  are  the  gain  matrix  obtained  from  the  Kalman 
filter  from  the  LQ  regulator. 

A  simulation  of  the  model  with  a  unit  impulse 
input  and  no  control  is  first  examined.  This  provides  a 
reference  to  refer  to  when  a  controller  is  later  added  to  the 
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system.  The  simulation  program  begins  with  an  interactive 
phase.  This  allows  the  user  to  select  the  values  of  certain 
parameters  that  are  not  fixed  by  either  the  program  or  input 
data.  The  parameters  that  may  be  varied  are: 

-  node  location  of  applied  disturbance 

-  axial  direction  of  the  disturbance 

-  sampling  time 

-  simulation  time 

-  damping  factor 

The  model  was  simulated  for  200  seconds  and  a  sampling  time  of 
one  second  which  was  based  on  the  period  of  the  highest 
sampling  frequency  of  the  system.  The  sampling  time  was 
approximately  ten  times  faster  than  the  highest  natural 
frequency,  resulting  in  a  minimal  amount  of  aliasing.  A 
damping  factor  of  0.1  was  used  in  order  to  yield  a  lightly 
damped  structure  and  mode  1  in  the  horizontal  translational 
direction  was  analyzed.  The  five  parameters  mentioned  above 
were  kept  constant  throughout  the  simulation. 

The  transient  response  of  both  mode  1  and  the 
total  displacement  due  to  an  impulse  input  is  shown  in  Figure 
10.  The  system  settles  in  approximately  800  seconds  when 
simulated  without  control.  The  ten  natural  modes  are  nearly 
equal  in  frequency.  As  a  result,  a  combination  in  which  2 
modes  reinforce  each  other  and  cancellation  may  then  result. 
This  is  called  the  beat  phenomenon. 
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Figure  10  Modal  Position  and  Displacement  with 
Control 


No 


Meirovitch  [Ref.  20]  states  that  any  motion  of  the  system  can 
be  regarded  at  any  time  as  a  superposition  of  the  natural 
modes  each  multiplied  by  some  constant.  The  contribution  of 
each  mode  is  represented  by  the  input  coupling  of  each 
vibrational  mode  as  shown  in  Figure  11.  It  can  be  seen  that 
modes  2  through  4  do  not  significantly  contribute  to  the 
response  of  the  system. 
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Figure  11  Input  Coupling  for  All  Modes 

The  total  energy  and  energy  per  mode  dissipated  by  the  system 
is  shown  in  Figure  12.  Although  the  system  response  deviates 
from  the  equilibrium  position,  the  total  energy  dissipated  by 
the  structure  is  minimal.  As  expected,  modes  2  through  4  have 
virtually  no  energy.  This  is  because  they  are  uncoupled  to 
the  input. 
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Figure  12  Total  Energy  and  Energy  Per  Mode  Dissipated 
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The  second  simulation  was  performed  with  the  LQ 
regulator  control  applied  to  the  system.  The  design  of  the  LQ 
regulator  for  the  LQG  controller  utilized  the  performance 
measure  in  Equation  3.5  with  R  equal  to  1.  The  stability  of 
the  system  is  crucial  to  the  evaluation  of  the  system 
performance  and  was  determined  by  evaluating  the  number  of 
encirclements  of  the  point  (0,-1)  which  represents  the  number 
of  poles  with  positive  real  parts.  This  evaluation  is  called 
the  Nyquist  criterion.  The  Nyquist  criterion  is  useful 
because  it  directly  displays  how  a  change  in  gain  affects 
stability  as  depicted  in  Figure  13.  The  second  plot  in  Figure 
13  is  an  enlargement  of  the  real  and  imaginary  axes  from  -2  to 
2.  Table  3  lists  the  matrix  gain  values,  K. 
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Table  3  Gain  Values 
Mode  1 


-7.8672E-02 
-8.0024E- 12 
-1 .0611E-11 
-2.3457E- 12 
4.8464E-02 
7.2167E-03 
-2.8356E-03 
6.6404E-02 
-1.3024E-03 


1 . 2409E+00 
-6.3944E- 1 1 
3 . 9083E - 1 0 
5.9271E-11 
-7.2182E-01 
-6. 1610E-01 
-2. 1306E-01 
8.1423E-01 
-1.4196E-02 


The  LQ  regulator  damps  the  oscillatory  motion  with  an  impulse 
input  as  illustrated  in  Figure  14.  The  transient  response 
settles  in  approximately  140  seconds. 
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Figure  14  Modal  Position  with  LQ  Regulator  for  Mode  1 


Figure  15  illustrates  the  total  energy  and  energy  per  mode  of 
the  system  controlled  by  the  LQ  regulator.  Using  full  state 
feedback,  the  energy  is  dissipated  rapidly  compared  to  Figure 
12  (without  control).  The  system  suppresses  the  vibratory 
motion  rapidly  and  the  higher  frequency  modes  damp  out  faster 
with  very  little  energy  lost. 


31 


Modes 


Figure  15  Total  Energy  Dissipated  with  LQ  Regulator 


Figure  16  illustrates  the  amount  of  control  energy  used  to 
suppress  the  oscillatory  motion.  Most  of  the  initial  control 
energy  was  used  for  mode  1  and  the  residual  for  higher 
frequency  modes.  It  was  observed  during  simulation  that  R  was 
inversely  proportional  to  the  control  energy  and  directly 
proportional  to  the  settling  time. 
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xtO'3  Control  Energy  fc  LQ  Regulator 


Figure  16  Control  Energy  for  LQ  Regulator 

The  LQ  regulator  provides  relatively  good  results  but  requires 
perfect  measurements  of  the  state.  Noise  is  inherent  in  any 
system  and  the  Kalman  filter  will  be  required  to  estimate 
these  states. 

b.  Simulation  of  Kalman  Filter 

The  addition  of  the  Kalman  filter  completes  the 
formulation  of  the  LQG  controller  as  depicted  in  Figure  4.  The 
parameters  to  be  determined  for  the  Kalman  filter  are  the 
scalar  plant  W  and  measurement  V  covariance  matrices.  This  was 
a  difficult  process.  The  sensor  or  measurement  noise 
covariance  matrix  value  was  determined  by  varying  its  value 
and  examining  the  system  response.  Several  values  for  W  were 
attempted  with  the  model.  Initially,  W=8002/3/40  and  V=0.005 
were  used.  The  value  for  W  was  too  high  and  required 
adjustment.  This  indicated  that  the  Kalman  filter  was 
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converging  much  faster  than  was  reasonable.  Consequently,  the 
measurement  noise  was  increased  and  W  was  lowered  which 
yielded  better  results.  It  was  determined  that  a  value  for  V 
=  0.005  and  W  =  40  provided  the  best  response  for  this  system 
as  show'n  in  Figure  17. 


Time  (sec) 

Figure  17  Modal  Position  and  Total  Displacement  with 
Kalman  Filter 


Figure  18  depicts  the  control  input  with  the  Kalman  filter. 
A  comparable  control  input  to  that  used  in  the  LQ  regulator 
was  required  to  suppress  the  vibrations. 
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Figure  18  Control  Input  with  Kalman  Filter 


Figure  19  illustrates  that  the  energy  dissipated  by  the  Kalman 
filter  is  essentially  the  same  as  that  of  the  LQ  regulator 
(Figure  15) . 


35 


o 

C 

Lu 


^  x1Q~*  ^  Toto'  Actuo!(  —  )ond  Estimatedf— .)£nengy  w/  K'olmqn  Filter 


0  20  40  60  80  100  120  140  160  180  200 


Time  (sec) 


Total  Energy  Per  Modg  with  Kalman  Filter 


~  1  ■ 

i . 

:  ■  1 

t  . 

i  i  : 

1 

. i . : . 1 

i 

j 

i  :  i  r— , _ 

.....  ,  i 

i 

r— i  i 

_ i _ i _ _ _ _ i 

4  e  e  io  i: 

Mo  dec 


Figure  19  Energy  Dissipated  by  Kalman  Filter 


Less  control  energy  was  required  early  in  the  response  with 
the  Kalman  filter  to  control  the  vibratory  motion  since  the 
estimates  have  not  yet  converged  to  the  states.  Figure  20 
shows  the  control  energy  required  by  the  LQG  controller. 
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X 1  :-4  Control  Enera»  for  LQ3  Controller 
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Figure  2  0  Control  Energy  with  LQG  Controller 

The  LQG  controller  not  only  guarantees  the  stability  of 
the  closed  loop  system,  but  also  minimizes  the  performance 
measure.  With  the  Kalman  filter  and  plant  input  noise  added, 
an  approximation  of  the  states  yielded  suitable  performance. 
Although  the  performance  of  the  system  is  relatively  good,  it 
may  no  longer  be  robust  to  perturbations  in  the  plant. 

To  increase  the  robustness  of  the  system,  loop 
transfer  recovery  was  applied.  As  a  result,  the  phase  and 
gain  margins  were  improved  to  the  values  indicated  in  Figure 
21.  The  second  plot  in  Figure  21  is  an  enlargement  of  the  real 
and  imaginary  axes  from  -2  to  2 .  A  compromise  was  made  in 
terms  of  the  robustness  and  performance  of  the  system. 
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Figure  21  Nyquist  Plot  for  Plant  with  Kalman  Filter 

Although  the  LQG  controller  is  the  optimal  controller  for  the 
modelled  noise  disturbances,  the  H®  design  provides  a 
controller  that  attenuates  disturbances  in  a  specified 
frequency  range  and  can  address  robustness  directly. 

4.  H«o  controller 

The  design  of  the  H®  controller  is  based  on  the 
desired  frequency  response  characteristics.  Dependent  upon 
this  configuration  is  the  selection  of  performance  and 
robustness  specifications  by  the  designer.  The  determination 
of  the  specifications  required  yields  the  necessary  weighting 
functions.  As  discussed  in  section  II. B. 4,  the  H®  controller 
design  is  analyzed  in  four  steps.  As  a  reminder,  the  four 
steps  are  to  formulate  a  problem  with  a  stable  feedback 
compensator,  determine  the  sensitivity  of  the  system,  augment 
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the  plant  with  weighting  functions,  and  apply  H=o  controller 
design  to  the  system. 

The  design  of  the  H°°  controller  is  based  on  the  space 
structure  example  in  Reference  21.  In  selecting  the  weighting 
functions,  the  0  dB  crossover  frequency  of  W.,(s)  must  be 
sufficiently  below  the  0  dB  crossover  frequency  of  W3(s).  This 
is  to  ensure  that  a  solution  exists.  W1  is  selected  as  a 
second-order  system  described  by: 


h\(s) 


y  ICO  (1  +  — —  ) : 
'  500 


(1  + 


(3.10) 


where  y  is  a  design  parameter.  W^s)  is  the  sensitivity 
specification  where  an  attenuation  of  100:1  to  10  rad/sec  (16 
hertz)  is  required.  The  robustness  specification  has  a  closed 
loop  bandwidth  100  rad/sec  (300  hertz).  W3(s)  is  selected  as 
a  derivative  function  of  the  first  order: 


W3  (s) 


s 

200 


(3.11) 


Figure  22  illustrates  that  the  requirement  of  the  crossover 
frequencies  is  met. 
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M*M3  Design  Specifications  —  Mode  7 
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Figure  22  Robustness  and  Sensitivity  Specification 

The  next  step  in  this  design  process  was  to  augment 
the  plant  matrix  to  generate  Figure  7.  The  maximum  singular 
values  of  the  sensitivity  and  complementary  sensitivity 
functions  must  be  less  than  one  where  the  closed  loop  system 
is  of  the  form: 


o  ( j  w) 


W.  ( j  w  ) 
W3  (jw) 


s  (jw) 
T(jw)  ) 


<;  l. 


(3.11) 


Once  the  weighting  functions  are  added  to  the  system,  the  H» 
controller  generated  acp,  bep,  ccp,  and  dep  in  state  space 
form: 

*C  =  acp  xc  +  bep  y 

uc  =  ccp  xc  +  dep  y  '  *  1 
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where  xc  is  the  controller  state  vector,  y  is  the  output,  and 
uc  is  the  control  input.  Combining  the  controller  equations 
and  the  plant  equations  yields  the  closed  loop  system  in 
Equation  3.13: 


’  x' 

A  +  C  B  dcp 

B  ccp 

X 

B 

*e. 

bcp  C 

acp 

Xc. 

0. 

The  inverse  of  both  weighting  functions  are  the  bounds  on  the 
frequency  response  of  the  sensitivity  and  complementary 
sensitivity  functions  of  the  closed  loop  system.  A  plot  of 
both  is  provided  in  Figures  23  and  24.  Gamma  equal  to  1.5  was 
used.  The  requirement  to  limit  the  sensitivity  to  not  exceed 
the  W1 " 1  was  met. 


1/*1  Sc  Sensitivity  Function  f©r  H-infinity  Controller 


Frequency  -  Rod /Sec 

Figure  23  1/Wl  &  Sensitivity  Function  for  H°° 
Controller 
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Frequenc*  -  R->:  Sr: 


Figure  24  1/W3  &  Complementary  Sensitivity  Function 
for  H°o  Controller 


The  controller  acts  as  a  filter  that  generates  the  control 
from  the  plant  output.  Applying  the  H°o  controller  resulted  in 
a  very  fast  closed  loop  system  as  shown  in  Figure  25: 
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Figure  25  Energy  and  Control  Input  with  H°o  Controller 


The  H°o  controller  designed  with  weighting  functions  W,(s)  and 
W3(s)  had  a  large  bandwidth  of  500  rad/sec  which  added  too 
much  control  to  the  system.  An  additional  weighting  factor 
was  added  to  the  performance  measure  to  lessen  the  control 
applied  to  the  system.  The  transfer  function  R(jw)  was 
constrained:  R(jw)  <  W2(jw).  R  has  no  common  name,  but  is 
defined  as: 


Uc(s) 

U^s) 


R(s)  =  F(  s)  [I  +  G(s)  K(s)  ]  ~ 1 


(3.14) 


where  R(jw)  is  the  ratio  of  the  control  input  to  the 
disturbance  input  in  the  real  frequency  domain.  W2(s)  was 
selected  to  be: 
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(3.15) 


10 9  (s  +  0.22) 
s  +  2 . 4el03 


There  was  no  significant  difference  as  a  result  of  this 
modification . 

The  next  modification  to  the  specifications  was 
lowering  the  bandwidth  of  the  system.  Several  values  for  the 
bandwidth  from  100  rad/sec  to  1  rad/sec  were  tried.  The 
frequency  range  from  1.25  rad/sec  to  2  rad/sec  provided  good 
response.  A  bandwidth  of  2,  however,  yielded  the  best  response 
and  is  shown  in  Figure  26. 
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Figure  26  Energy  and  Control  Input  with  H®  Controller 

Figure  26  illustrates  that  less  control  energy  and  control 
input  was  required  to  actively  suppress  the  vibrations. 
Moreover,  mode  1  required  the  most  control  and  was  suppressed 
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immediately  with  residual  control  for  the  other  modes.  The 
total  energy  in  the  system  dissipates  in  about  20  seconds. 

The  full  order  design  for  the  H®  controller  was  based 
on  trial  and  error.  A  better  design  may  exist.  The  full  order 
controller  yielded  good  results.  The  total  energy  to  suppress 
the  vibrations  with  the  H®  controller  was  greater  than  that 
with  LQG  controller.  However,  it  required  only  20  seconds  vice 
140  seconds.  The  full  order  controllers  require  more  CPU  time, 
hence  reduced  order  controllers  are  implemented  and  examined. 

B.  REDUCED  ORDER  CONTROLLERS 

1.  Linear  Quadratic  Gaussian  Controller 

The  reduced  order  models  are  obtained  by  truncating 
the  high  frequency  modes.  The  performance  measure  and  design 
parameters  remain  the  same  as  those  for  the  full  order 
controller.  Controllers  are  designed  for  the  reduced  order 
models  and  applied  to  the  full  order  system.  The  simulation 
software  allows  the  user  to  specify  the  number  of  modes  to  be 
controlled  through  an  interactive  phase  where  the  user  is 
queried  for  the  number  of  modes  to  be  controlled.  As  the 
number  of  modes  to  be  controlled  decreases,  less  energy  is 
dissipated  as  shown  in  Figure  27. 
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Figure  27  Total  Energy  Dissipated  with  7  Modes 
Controlled 


Figure  28  shows  that  the  control  energy  also  decreased.  This 
is  because  there  are  fewer  modes  contributing  to  the  control. 
The  more  modes  controlled  the  closer  the  reduced-order 
controller  was  to  the  full-order  design. 
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Figure  28  Control  Energy  with  7  Modes  Controlled 

With  the  Kalman  filter  added,  the  LQG  controller 
exhibits  good  response  with  the  impulse  input.  The  actual  and 
estimated  energy  with  the  reduced-order  Kalman  filter  is  shown 
in  Figure  29: 
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Figure  29  Reduced  Order  Controller  with  7  Modes 
Controlled 

The  estimation  of  the  actual  position  required  approximately 
the  same  amount  of  time  as  Figure  18  (full  order  controller). 

2.  H«  controller 

The  design  of  the  H°°  order  reduced  controller  is 
similar  to  that  of  the  LQG  controller.  A  bandwidth  of  2 
rad/sec  was  used  for  the  system  again.  Figure  30  illustrates 
damping  of  the  modes  and  shows  the  amount  of  control  and 
control  energy  with  7  modes  controlled: 
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Figure  30  Energy  and  Control  Input  for  H°°  Reduced 
Controller 

The  performance  of  the  reduced-order  system  was  reasonable  and 
used  less  energy  than  the  full  order  controller. 

C.  A  COMPARISON  OF  LQG  AND  H«  CONTROLLERS 

Both  design  processes  involved  some  insight  on  behalf  of 
the  designer  to  formulate  the  performance  measures,  i.e.,  the 
plant  and  measurement  noise  covariance  matrices  for  the  LQG 
controller  and  the  weighting  functions  for  the  H®  controller. 
The  LQG  controller  exhibited  good  stability  margins  and  also 
minimized  the  performance  measure  specifying  by  two  weighting 
functions  Q  and  R.  The  H®  controller  specified  the 
performance  measure  using  the  weighting  functions,  W^s)  and 
W3(s).  Thus,  both  controllers  required  the  selection  of 
weighting  functions. 
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As  expected  the  L QG  controller  provided  good  full  order 
response  with  modelled  noise.  For  robustness,  another 
controller  may  be  better  suited  because  no  guidelines  exist  to 
tailor  the  LQG  controller  for  unmodelled  dynamics.  The  H°° 
controller  is  designed  for  such  a  condition. 

The  full-order  controller,  of  course,  exhibited  better 
response  than  the  reduced  order  controller.  However,  the 
reduced-order  controllers  worked  reasonably  and  used  less 
energy  for  control. 

An  indication  of  the  robustness  of  the  LQG  and  H°o 
controller  is  the  performance  of  the  reduced-order 
controllers.  Both  controllers  proved  to  be  reasonably  robust 
to  unmodelled  dynamics. 
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IV.  CONCLUSIONS 


The  significant  observations  and  results  from  the 
simulations  generated  are  summarized.  Areas  that  require 
further  consideration  for  research  and  development  are 
highlighted . 

A.  SUMMARY  OF  OBSERVATIONS 

The  design  of  the  full  order  controllers  was  based  on  the 
full  order  model  of  the  system.  As  a  result,  the  "best" 
performance  of  the  system  was  expected.  The  LQG  controller 
illustrated  good  performance  in  the  presence  of  modelled  noise 
dynamics.  However,  it  lost  some  of  its  robustness  properties 
due  to  estimxting  the  states  with  the  Kalman  filter.  Applying 
LTR,  the  robustness  was  increased  by  implanting  fictitious 
noise,  thereby  yielding  a  robust  suboptimal  system.  The  H°o 
controller  is  an  optimal  method  directly  designed  for 
robustness  to  unmodelled  dynamics. 

The  H°°  controller  design  is  simple  and  yet  complicated  for 
implementation.  It,  like  the  design  of  the  LQG  controller, 
requires  some  perspective  in  its  design.  The  weighting 
functions  used  for  augmentation,  provide  a  tool  for  specifying 
robustness.  The  weighting  functions  characterize  the 
specifications  on  disturbance  attenuations  and  robustness. 
Simulations  of  the  system  with  both  the  LQG  and  H°o  controller 
were  then  performed. 
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Simulation  results  indicated  that  both  dynamic  controllers 
provided  good  results.  The  LQG  controller  provided  good 
response  with  less  control  energy  required  than  the  Hoc 
controller.  The  reduced  order  LQG  and  Hoo  controllers  both 
achieved  good  response  when  truncating  the  high  frequency 
modes.  Simulations  of  the  system  with  7  modes  controlled  were 
presented.  It  was  shown  that  the  reduced  order  controller 
actively  suppressed  the  modes  selected  to  be  controlled.  The 
reduced  order  controllers  used  less  energy  and  control  as 
expected . 

B.  FURTHER  RESEARCH  AND  DEVELOPMENT 

1.  The  range  of  natural  frequency  values  should  be 
increased  in  order  to  better  evaluate  the  system  response. 

2.  The  model  should  be  evaluated  with  random  noise  input 
disturbances . 

3.  /i-synthesis  should  be  applied  to  the  model  and  then 
compared  to  the  LQG  and  H®  controllers. 
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Appendix  A 

Zeros  and  Poles  of  Plant 


Zeros 

-4.5920e  +  07-  1.2343e  +  05  i 
4.5920e  +  07+1.2343e  +  05i 
.2193e-03-4.4395e-01i 
.  1 600e-03-4 .3 1 73e-0 1  i 
.0590e-03-4. 1223e-0 1  i 
.963  le-03-3.9319e-01i 
.0148e-03-4.0345e-01i 
.2200e-03  +  4.4407e-01i 
.  1 600e-03  +  4.3 1 73e-0 1  i 
.0590e-03  +  4.1223e-01i 
.0148e-03  +  4.0345e-01  i 
.0391e-03  +  4.0753e-01  i 
-1.963  le-03  +  3.93 19e-01i 
-1.7350e-03-3.4727e-01i 
- 1 .8900e-03-3.788 1  e-0 1  i 
-1.9050e-03-3.8105e-01i 
-1.7350e-03  +  3.4727e-01i 
-1.8900e-03  +  3.788  le-Oli 
- 1 .9050e-03  +  3.8 105e-01  i 
-2.21 93e-03  +  4.4395e-0 1  i 


Poles 

-1.6200e-03  +  3.2372e-01 
-1.6200e-03-3.2372e-01 
-1.7350e-03  +  3.4727e-01i 
-1.7350e-03-3.4727e-0li 
-1.8900e-03  +  3.7881e-01i 
-1.8900e-03-3.788  le-Oli 
-1.9050e-03  +  3.8105e-01i 
- 1 .9050e-03-3.8 1 05e-0 1  i 
-1.9800e-03  +  3.9619e-01i 
-1.9800e-03-3.9619e-0li 
-2.0250e-03  +  4.0533e-0 1  i 
-2.0250e-03-4.0533e-0 1  i 
-2.0400e-03  +  4.0767e-0 1  i 
-2.0400e-03-4.0767e-01i 
-2.0750e-03  + 4. 152  le-Oli 
-2.0750e-03-4. 152  le-Oli 
-2. 1 600e-03  +  4.3 1 74e-0 1  i 
-2. 1 600e-03-4.3 1 74e-0 1  i 
-2.2200e-03  +  4.4407e-0 1  i 
-2.2200e-03-4.4407e-01  i 
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Appendix  B 
Computer  Programs 


************************************************************ 

*  LQGCTC.M 

*  This  program  generates  the  Linear  Quadratic  Regulator  for  the 

*  the  system.  * 

(•ox****************************************************** 


%  In  order  to  run  this  program  some  of  the  variables  used  in  the 
9 c  model  simulation  program  are  required. 

if  exist('modenbr')  =  =  1 
load  lqgcc.mat 
else 
var 
end 

9cFor  Linear  Gaussian  Control,  the  system  is  modeled  as: 

%  aida_dot  =  A*aida  +  Blu  +  B2w 
9  c  y  =  C*aida  +  v 

c/c  where  w  is  the  plant  driving  noise  and  v  is  the  measurement  noise 
9 c  Control  input  is  u=  -K*aida  where  aida  is  the  modal  amplitude. 

9/t  Using  the  ten  second  order  matrices,  the  response  should 
9r  be  quicker  than  that  of  the  twentieth  order . 

%  Plant  matrices  for  20th  order  system  to  obtain  full  order 
9cstate  feedback  controller 

dv  =  zeros(19, 1); 
for  i  =  1:19 
if(rem(i,2)'  =0); 

dv(i,l)=  1; 
end 
end 


A  =  diag(dv,  1 )  +  diag([-0. 1  W8,0,-0. 1206,0,-0. 1435,0,-0. 1452,0,... 

-0. 1 5697,0,-0. 1643,0,-0. 1 662,0,-0. 1 724,0,-0. 1864,0,-0. 1 972],- 1 )  + ... 
diag([0, -0.00324,0, -0.00347, 0,-0.00378,0,-0.00381,0,-0.00396,0,... 

-0.00405,0,-0.00408.0,-0.004 1 5,0, -0.00432, 0,-0.00444],0); 

B  =  [0;b(l);0;b(2);0;b(3);0;b(4);0;b(5);0;b(6);0;b(7);0;b(8);0;b(9);... 
0;b(  10)]; 

Q  =  diag([f(l),l,f(2),l,f(3),l,f(4).l,f(5),l,f(6),l,f(7),l,f(8),... 
l,f(9),l,f(10),l],0); 


R=  l; 


9c  Full  Order  Stale  Feedback  where  u  =  -K*[aidal,aida2,aida3,...]' 

9c  Phi  represents  the  discrete  2x2  matrix  for  each  mode  of  vibration 
9c  at  natural  frequency.  These  matrices  can  be  found  on  the  diagonals 
9c  of  this  twentieth  order  matrix.  Del  also  represents  each  mode  of 
%  vibration.  First  two  row  entries  correspond  to  a  modal  matrix 
9c  at  specified  natural  frequency. 

ts  =  1 ;  9c Sampling  Time 

[Phi,Del]  =  c2d(A,B,ts); 

aida(:,l)  =  B;  9 c  Initialization 

K  =  dlqr(Phi,Del,Q,R); 

for  n  =  1:200  Cc  Time 

u  =  -K*aida(:,n);  9c  No  control  added 

for  mode  =1:10 

Phimode  =  Phi(2*  mode- 1:2*  mode, 2*  mode- l:2*mode); 

Deimode  =  Del(2*mode-l:2*mode); 

Qe  =  Q(2*mode-l:2*mode,2*mode-l:2*mode);  9c  Energy  Q  matrix 

9c  Aida  gives  the  matrix  for  twenty  states  for  each  point  in  time. 

9c  Next  actual  state  calculated 

aida(2*mode-l:2*mode,n  +  1)  =  Phimode*aida(2*mode-l:2*mode,n)... 

+  Deimode  *  u  ;  %  control  added 

9c  At  each  point  in  time  the  energy  of  each  mode  is  calculated  using 
9c  Energy  in  specified  mode  =  0.5[aida  aida_d]0[aida 
9c  aida_d] 

E(mode,n)  =  0.5*[aida(2*mode-l,n)  aida(2*mode,n)]  *... 

Qe  *[aida(2* mode- l,n);aida(2* mode, n)]; 

end 

end 

%  Plot  of  Total  Dissipation  of  Modes 
for  m  =  1:10 

out_d(:,m)  =  aida(2*m-l,:)’*B(2*m); 
end 

outtot  =  sum(out_d'); 

9c  Plot  of  Total  Displacement 
time  0  =  0:1:200; 
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plot(time_o,out_tot); 

xlabel(Time(sec)  ');ylabel('Amplitude');grid; 
title(Total  Displacement  of  Modes'); 
gtext(['R  =',  num2str(R)  ]); 
clg; 

%  Energy'  Calculation 
E_total  =  sum(E); 
time  =  0:l:199; 

subplot(211),plot(E_total);grid; 
title(Total  Energy  with  Control’); 
xlabel(’Time'); 

gtext(['R  =  ',num2str(R)]);pause 

%  Plot  of  Total  Energy'  per  Mode 

Em_total  =  sum(E’); 

modes=  1:1:10; 

plot(modes,Em_total);grid; 

title(’Total  Energy  per  Mode  with  Control'); 

gtext(['R  =  ',num2str(R)]); 
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************************************************************************ 


*  KALMAN.  M 

*  This  program  generates  the  Kalman  filter  for  the  LQG  controller.  * 

************************************************************************ 


c/c  In  order  to  run  this  program  some  of  the  variables  used  in  the 
%  model  simulation  program  are  required.  This  program  simulates  the  model 
9c  with  estimated  states  using  steady  state  Kalman  filter. 


clg 

format  short  e 
format  compact 

if  exist('b')  =  =  1 
load  var.mat 
else 
var 
end 

%  Plant  matrices 
dv  =  zeros(19,l); 
for  i  =  1:19 
if(rem(i,2)'  =0); 

dv(i,l)=l; 

end 

end 


A  =  diag(dv.  1 )  +  diag([-0. 1048,0,-0. 1206,0,-0. 1435,0,-0. 1452,0,... 

-0. 15697,0,-0. 1643,0,-0. 1 662,0,-0. 1724,0,-0. 1 864,0,-0. 1 972],- 1 )  + ... 
diag([0,-0.00324,0,-0.00347,0,-0.00378,0,-0.0038 1,0,-0.00396,0,... 

-0.00405,0,-0.00408,0, -0.00415, 0.-0.00432,0, -0.00444], 0); 

B  =  [0;b(l);0;b(2);0;b(3);0;b(4);0;b(5);0;b(6);0;b(7);0;b(8);0;b(9);... 

0;b(  10)]; 

C  =  [b(  1)  0  b(2)  0  b(3)  0  b(4)  0  b(5)  0  b(6)  0  b(7)  0  ... 

b(8)  0  b(9)  0  b(  10)  0]; 

D  =  0; 

O  =  diagf  [f(  1 ),  1  ,f(2),  1  ,f(3),  l,f(4),  1  ,f(5),  l,f(6),  l,f(7),  l,f(8),... 
l,f(9),  l,f(  10),  1  ],0); 

R=  l; 

N  =  200; 

9c  Initialize  the  random  inputs  to  the  same  for  each  run. 
rand('normal’); 

randfseed',0);  9c  Sets  the  seed  to  0  when  Matlab  is  entered 
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%  Initialization  of  the  state  and  the  estimate,  x_hat 
%  and  the  error  covariance  matrix. 
aida_hat(:,l)  =  zeros(B); 
aida(:,l)  =  zeros(B); 
old_P  =  zeros(A); 

%  Desired  state 
aida_d  =  zeros(B); 


ts  =  1;  %Sampling  Time 

W  =  input(Tnput  W:  '); 

V  =  input(’Input  V:  ’); 

[Phi, Del]  =  c2d(  A,B,ts); 

G  =  real(dlqe(Phi,Del,C,W,V));  %  Steady  State  Optimal  Gain 

K  =  real(dlqr(Phi,Del,Q,R)); 

for  k=  1:N 

c/c  Plant  simulation 
ww  =  zeros(N,l);ww(l)  =  1; 
u(k)  =  -K*(aida_hat(:,k)-aida_d); 
aida(:,k+  l)  =  Phi*aida(:,k)  +  Del*u(k)  +  Del*ww(k); 
w(k+  l)  =  sqrt(V)*rand; 
y(k+ l)  =  C*aida(:,k+ 1)  +  w(k+l); 

%  Estimates  updated 

aida_hat(:,k+ 1)  =  Phi*aida_hat(:,k)  +  Del  *u(k); 
y_hat(k+ 1)  =  C  *aida_hat(:,k+  1); 

aida_hat(:,k+ 1)  =  aida_hat(:,k+  1)  +  G*(y(k  + l)-y_hat(k+ 1)); 
t(k)  =  k-l; 


for  mode  =  1:10 

E_k(mode,k)  =  0.5 *[aida(2* mode- l,k)  aida(2*mode,k)]*... 
Q(2*mode-l:2*mode,2*mode-l:2*mode)*... 

[aida(2*  mode- l,k);aida(2*  mode, k)]; 
end 
end 

u(N+  l)  =  -K*aida_hat(:,N+  1); 
t(N+  1)  =  N+  1; 

%  Plots 
subplot(21 1) 

plot(t,aida(modenbr,:),’*',t,aida_hat(modenbr,:),'o’); 
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title([’Estimated  and  Actual  Position  of  Mode  \num2str(modenbr)]); 

xlabelCseconds'^ylabelCDisplacement'^grid; 

pause 

subplot(212),plot(t,u);grid; 
title('Control  Input');pause 
%meta  kalman 
clg 

%  Plot  of  Total  Dissipation  of  Modes 
for  d  =  1:10 

k_disp(:,d)  =  aida(2*d-l,:)’*B(2*d): 
kh_disp(:,d)  =  aida_hat(2*d-l.:)'*B(2*d); 
end 

t  disp  =  sum(k_disp’); 
th_disp  =  sum(kh_disp'); 

%  Plot  of  Total  Displacement 
ktime  =  0:l:N; 

plot(ktime,t_disp,'--',ktime,th_disp, 

xlabel('Time(sec)  '):ylabel('Amplitude');grid; 

title('Total  Displacement  of  Modes  with  Kalman  Filter  ’); 

%meta  kalman 
pause 

%  Plot  of  Total  Energy 
K_energy  =  0.5*(aida'*Q*aida); 

Khenergy  =  0.5*(aida_hat'*Q*aida_hat);clg; 

plot(ktime,diag(K_energy),'--',ktime,diag(Khenergy),'-.’);grid; 

title(Total  Actual(--)and  Estimated(-.)Energy  w/  Kalman  Filter’); 

xlabel(Time');ylabel('Magnitude'); 

pause; 

%meta  kalman 
clg 

%  Plot  of  Total  Energy  per  Mode 
modes  =  1:1:10; 

Emk_tot  =  sum(E_k');clg; 
plot(modes,Emk_tot);grid; 

title(Total  Energy  Per  Mode  with  Kalman  Fiiter');xlabel('Modes  '); 

ylabel('Amplitude');pause 

%meta  kalman 


59 


************************************************************************ 


»  JOSE.M 

*  This  program  generates  the  H»  controller  for  the  system.  This  program* 

*  has  been  modified  to  reflect  the  specification  design  for  * 

*  this  particular  model.  * 

************************************************************************ 


%  JOSE  Large  space  structure  design  demonstration 

% 

%  R.  Y.  Chiang  &  M.  G.  Safonov  3/88 
%  Copyright  (c)  1988  by  the  MathWorks,  Inc. 


%  All  Rights  Reserved. 

% . 

clc 

disp('  <  <  Demo  #3:  MIMO  Large  Space  Structure  Design  Example  >  >') 

dispC  _ ') 

disp('  Secondary  Mirror -—>  -ooo-  A’) 

dispC  /\/\  D 

dispC  I  X  |  |') 

disp('  I  /  \  i  D 

disp('  .  i') 

disp('  |  \  /  |  I') 

dispC  |  X  |  I’) 

dispC  i  /  \  I  ’) 

disp(’  Lens . >  — O—  7.4  Meters') 

disp('  |  |  ’) 

dispC  I  \  /  I  D 

dispC  I  \  /  I  I') 

disp('  I  X  |  n 

dispC  |  /  \  |  D 

dispC  I  /  \  I  I') 

dispC  II  \  \  I') 

dispC  Primary  Mirror  ->  (  OOOOOO _ \ _ )  |’) 

disp('  \  /I') 

dispC  \ _ / _ v _ ') 

dispC  ') 

dispC  (strike  a  key  to  continue  ...)') 

pause 

if  exist  (’b')==  1 


load  var.mat 

lse 

var 


60 


format  short 
format  compact 


dv  =  zeros(19,l); 
for  i  =  1:19 

if  (rem(i,2)*  =0); 

dv(i,l)  =  1; 
end 
end 

ag  =  diag(dv,l )  +  diag([-0. 1048,0,-0. 1206,0,-0. 1435,0,-0.1452.0,... 

-0. 15697,0,-0. 1 643,0,-0. 1 662,0,-0. 1 724,0,-0. 1 864,0,-0. 1 972],- 1 )  + ... 
diag(  [0,-0.00324, 0,-0.00347,0,-0.00378, 0,-0.00381, 0,-0.00396, 0,... 

-0.00405,0,-0.00408.0, -0.00415, 0,-0.00432,0,-0.00444],0); 
bg  =  [0;h(l);0;h(2);0;b(3);0;b(4);0;b(5);0;b(6);0;b(7);0;b(8):0;b(9);... 
0;b(  10)]; 

cg  =  [b(l)  0  b(2)  0  b(3)  0  b(4)  0  b(5)  0  b(6)  0  b(7)  0  ... 

b(8)  0  b(9)  0  b(  10)  0]: 
dg  =  0; 
clc 

dispC  ') 

disp('  State-Space  of  the  Large  Space  Structure:') 

dispC  (after  colocated  rate  feedback  and  model  reduction)') 

ag 

bg 

dispC  (strike  a  key  to  continue  ...)’) 

pause 

clc 

eg 

dg 

disp('  (strike  a  key  to  continue  ...)') 

pause 

clc 

disp(’  ’) 
disp(’  ’) 

dispC  Poles  of  the  Plant  (stable):') 
dispC  ') 


dispC  . - . - — ’) 

disp('  poleg  =  eig(ag)  ck  Computing  the  poles  of  the  plant') 

dispC  . - . ') 

poleg  =  eig(ag) 


pause 
dispf  ') 
dispC  ') 
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(strike  a  key  to  continue  ...)') 


disp(' 
pause 
clc 

disp(’  ') 
dispC  ') 

dispC  Transmission  Zeros  of  the  Plant  (minimum  phase):') 


disp('  ') 

disp('  . . . ') 

dispC  tzerog  =  tzero(ag.bg,cg,dg)  9 c  Computing  the  transmission  zeros' 

dispC  . ') 

tzerog  =  tzero(ag,bg,cg,dg)  9 c  Computing  the  transmission  zeros 

dispC  ') 
dispC  ') 
dispC  ') 

dispC  (strike  a  key  to  continue  ...)') 

pause 

clc 

dispC  ') 

dispC  - Computing  SV  Bode  plot  of  the  open  loop  plant  -  -  -') 


w  =  logspace(-3,5,100); 

svg  =  sigma(ag,bg.cg,dg,l,\v);  svg  =  20*logl0(svg); 

dispC  ') 
dispC  ’) 
dispC  ') 
disp('  ') 
dispC  ') 

dispC  (strike  a  key  to  see  the  SV  Bode  plot  of  G(s)  ...)’) 

pause 

semilogx(w,svg) 

title(’MIMO  --  Mode  1:  Open  Loop  Bode  Plot  ’) 

xlabel('Frequency  -  Rad/Sec') 

ylabel('SV  -  db') 

grid 

pause 

%  meta  bode 
clc 

dispC  ') 

dispC  <  <  Design  Specifications  >  >  ') 

dispC  ’) 

dispC  1).  Robustness  Spec.  :  closed  loop  bandwidth  --  200  r/s  (30  hz)') 
disp('  Associated  Weighting:') 

dispC  ') 

dispC  -1  200') 

dispC  W3(s)  =  .  *  I  (fixd)') 
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2x2') 


dispC  s 

dispC  ') 
dispC  ') 

dispC  2).  Performance  Spec.:  sensitivity  reduction  of  at  least  100:1') 
dispC  up  to  approx.  100  r/s’) 

dispC  Associated  Weighting:') 

dispC  •) 

dispC  -1  -1  0.01(1  +  s/10)  “2') 

dispC  Wl(s)  =  Gam  *  - *  I’) 

dispC  (1  +  s/500) '2  2x2') 

dispC  ') 

dispC  where  "Gam"  in  our  design  goes  from  1  —  >  1.5.') 

coef=input('Enter  the  coef  value'): 
fac  =  500/coef; 
k  =  200/fac; 

nuw3i  =  [0  k];  dnw3i  =  [1  0]; 

svw3i  =  bode(nuw3i,dnw3i,w);  svw3i  =  20*logl0(svw3i); 

nuwli  =  conv([l/(10/fac)  l].[l/(  10/fac)  lj); 

dnwli  =  100*conv([l/coef  l],[l/coef  1]); 

svwli  =  bode(nuwli,dnwli,vv);  svwli  =  20*logl0(svwli); 

dispC  ') 

dispC  ') 

dispC  (strike  a  key  to  see  the  plot  of  the  weightings  ...)') 

pause 

axis([0  5  -40  40]) 
semilogx(w,svw  1  i,w,svw3i ) 
grid 

title('MlMO  Design  Specifications  -  Mode  1') 
xlabel('Frequencv  -  Rad/Sec') 
ylabel('l/Wl  &  1  /\V3  -  db') 
text(2, -20, ’Sensitivity  Spec.-  1/Wl(s)’) 
text(2, 10, 'Robustness  Spec.-  1/W3(s)’) 
pause 

%  meta  weight 

axis 

clc 

dispC  <<  Problem  Formulation  >>') 

dispC  ') 

dispC  Form  an  augmented  plant  P(s)  with  these  two  weighting  functions:’) 
dispC  ') 

dispC  1).  Wl  penalizing  error  signal  "e"') 

dispC  ’) 

dispC  2).  W3  penalizing  plant  output  "y"') 

dispC  ') 
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disp(' 

disp(’ 

dispC  ') 

disp(' 

disp(' 

disp('  ') 

disp(' 

dispC 

dispC 

disp(' 

dispC 

dispC  ') 

dispC  ') 

dispC  ') 

dispC 

pause 

clc 

dispC  ') 

dispC  ') 

dispC 

dispC  ') 

dispC 

disp(' 

dispC 

dispC 

dispC 

dispC 

dispC 

disp(' 

dispC 

dispC 

dispC 

dispC 

dispC  ') 

dispC  ’) 

dispC 

pause 

clc 

dispC  ’) 
dispC  ') 
dispC 
dispC 
dispC 
dispC  ’) 


and  find  a  stabilizing  controller  F(s)  such  that  the  Hinf-norm’) 
of  TF  Tylul  is  minimized  and  less  than  one,  i.e.') 


min  |Tylul  |  <  1,') 

F(s)  '  inf) 


where  ') 


i  -in 

Tvlul  =  |  Gam*Wl*(I  +  GF)  |  =  |  Gam 
i  -ii  i  W3  *  (i  - s)  n 

|  W3*GF*(I  +  GF)  j’) 


W1 


S  |') 


(strike  a  key  to  continue  ...)’) 


<  <  DESIGN  PROCEDURE  >  >  ) 


*****************************> 

[Step  1],  Do  plant  augmentation  (run  augtf.m  or 


augss.m) 


’) 


") 


[Step  2].  Balanced  the  augmented  plant  for  better  *') 

numerical  condition  (run  OBALREAL.M)  *') 

*') 

[Step  3].  Do  H-inf  synthesis  (run  HINF.M)  *') 

*') 

[Step  4],  Redo  the  plant  augmentation  and  balancing  *') 
for  a  new'  "Gam"  -->  1.5  and  rerun  HINF.M  *’) 

') 


****************************  **\ 


(strike  a  key  to  continue  ...)’) 


’) 


Assign  the  cost  coefficient  "Gam"  -->  1  ') 
this  will  serve  as  the  baseline  design  ....') 
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dispC  ') 

Gam  =  input('  Input  the  cost  coefficient  "Gam"  =  '); 

dispC  ') 


disp('  . ') 

dispC  %  Plant  augmentation  of  the  LSS:') 

dispC  sysg  =  [ag  bg;cg  dg];  xg  =  4;') 

dispC  wl  =  [Gum*dnwli;nu\vli;Gam*dnwli;nuwli];’) 

dispC  w2  =  [];’) 

dispC  w3  =  [dnw3i;nuw3i;dnw3i;nuw3i];’) 

dispC  [A,B1,B2,C1,C2,D11,D12,D21,D22]  =  augtf(sysg,xg,wl,\v2,w3);') 

dispC  . . ') 


sysg  =  [ag  bg;cg  dg]:  xg  =  20; 
wl  =  [Gam*dnwli;nuwli]; 

%nuw2i  =  10e9’[l  0.12];dnw2i  =  [1  2.4e3]; 

%w2  =  [dnw2i;nuw2i]; 
w2  =  []; 

w3  =  [dnw3i;nuw3i];  %dnw3i;nuw3i]: 

[A,Bl,B2,Cl.C2,Dll,D12,D21,D22]  =  augtf(sysg,xg,wl,w2,w3); 
dispC  ’) 


dispC  -  -  -  State-Space  (A,B1,B2,C1,C2,D11,D12,D21,D22)  is  ready  for') 

disp('  the  Small-Gain  problem - ') 

dispC  ') 

%disp('  . ') 

%disp('  [aa,bb,cc,mm,tt]  =  obalreal(A,[Bl  B2],[C1;C2])  %  Balancing  P(s)’) 
%disp(’  A  =  aa;  B1  =  bb(:,l:2);  B2  =  bb(:,3:4);  ') 

%disp('  Cl  =  cc(l:4,:);  C2  =  cc(5:6,:);') 

%disp('  . - . ') 


%[aa,bb,cc,mm,tt]  =  obalreal(A,[Bl  B2],[C1;C2]); 

%A  =  aa;  B1  =  bb(:,l:2);  B2  =  bb(:,3:4);  Cl  =  cc(l:4,:);  C2  =  cc(5:6,:); 
dispC  ') 
dispC  ’) 

dispC  (strike  a  key  to  continue  ...)') 

pause 

clc 

disp(’  ’) 
dispC  ') 
dispC  ') 


dispC  . - . - . ') 

dispC  hinf  %  Running  script  file  HINF.M  for  H-inf  optimization’) 

dispC  . ') 

hinf 

dispC  ') 
dispC  ') 

dispC  (strike  a  key  to  continue  ...)') 
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pause 

pltopt  %  Preparing  singular  values  for  plotting 

svwlil  =  svwli;  hsvsl  =  svs;  hsvtl  =  svt;  hsvttl  =  svtt; 
dispC  ') 
dispC  *) 

dispC  (strike  a  key  to  continue  ...)’) 

pause 


clc 

disp('  *) 
dispC  ’) 

dispC  After  a  few  iterations,  we  found  a  new  Gam  of  1.5  can  push  the') 
dispC  ') 

dispC  H-inf  cost  function  close  to  its  limit.  ') 
dispC  ') 
dispC  ') 

dispC  Input  "Gam"  -->  1.5,  and  try  HINF  again . ’) 

dispC  ') 
dispC  ’) 

Gam  =  input('  Input  the  cost  coefficient  "Gam"  =  '): 

dispC  ') 

dispC  . ') 

dispC  %  Adjust  plant  augmentation:') 

dispC  wl  =  [Gam*dnwli;nuwli;Gam*dnwli;nuwli];') 

dispC  [A,B1,B2,C1,C2,D11,D12,D21,D22]  =  augtf(sysg,xg,wl,w2,w3);’) 

dispC  . ') 

wl  =  [Gam*dnwli;nuwli];  %Gam*dnwli;nuwli]; 


[A,B1,B2,C1,C2,D1 1,D12,D21,D22]  =augtf(sysg,xg,wl,w2,w3); 
%[aa,bb,cc,mm,tt]  =  obalreal(A,[Bl  B2],[C1;C2]); 

%A  =  aa;  B1  =  bb(:,l:2);  B2  =  bb(:,3:4);  Cl  =  cc(l:4,:);  C2  =  cc(5:6,:); 
dispC  ') 
disp('  ') 

dispC  (strike  a  key  to  continue  ...)') 

pause 

hinf 

dispC  ') 
dispC  ') 

dispC  (strike  a  key  to  continue  ...)') 

pause 

pltopt 

svwli2  =  svwli;  hsvs2  =  svs;  hsvt2  =  svt;  hsvtt2  =  svtt; 
disp('  ’) 
dispC  ’) 

disp('  (strike  a  key  to  see  the  plots  of  the  comparison  ...)') 

pause 
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se  milogx(w,svw  1  i  1  ,w,  hsvs  1  ,w,svw  1  i2,w,hsvs2 ) 

title('Mode  1:  1/W1  &  Sensitivity  Function  for  H-infinity  Controller') 
xlabel('Frequencv  -  Rad/Sec’) 
ylabeK'SV  -  db')’ 
grid 

text(0.002,0,'H-inf  (Gam  =  1)  — >  H-inf  (Gam  =  1.5)') 

%  meta  sens 
pause 

semilogx(w,svw3i,w,hsvtl,w,hsvt2) 

title('  1/W3  &  Comp.  Sensitivity  Function  for  H-infinity  Controller') 
xlabel('Frequency  -  Rad/Sec') 
ylabel('SV  -  db')’ 
grid 

text(0.002,0, 'H-inf  (Gam  =  1)  — >  H-inf  (Gam  =  1.5)’) 

%  meta  compl 
pause 

semilogx(w,hsvttl,w,hsvtt2) 

title('Mode  1:  H-infinity  Design  Cost  Function  ') 

xlabel('Frequency  -  Rad/Sec') 

ylabel('SV  -  db') 

grid 

text(0.002, -10, 'H-inf  (Gam  =  1)  — >  H-inf  (Gam  =  1.5)') 

%  meta  cost 

pause 

clc 

disp('  ') 
disp('  ') 

disp('  <  <  8-State  H-inf  Controller  (Gam  =  1.5)  >  >') 

dispO  ') 

disp('  Poles  of  Controller  :') 
polecp  =  eig(acp) 
disp('  ') 

disp('  (strike  a  key  to  continue  ...)') 

pause 

clc 

disp('  State-Space  of  the  8-State  H-inf  Controller:') 
disp('  First  6  columns  of  the  A  matrix:') 
acp(:,l:6) 

disp('  (strike  a  key  to  continue  ...)') 

pause 

clc 

disp('  ’) 

disp('  Last  two  columns  of  the  A  matrix:') 
acp(:,7:8) 
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(strike  a  key  to  continue  ...)') 


dispC 
pause 
clc 
bcp 

disp(’  (strike  a  key  to  continue  ...)') 

pause 

clc 

ccp 

dcp 

disp('  (strike  a  key  to  continue  ...)') 

pause 

clc 

disp('  ’) 

disp('  Poles  of  closed-loop  TF  matrix  Tvlul:') 
poletyu  =  eig(acl) 

disp('  (strike  a  key  to  continue  ...)’) 

pause 

r-f 

VC 

save  hinfcc.mat  ag  bg  eg  dg  acp  bcp  ccp  dcp  b  f  modenbr  A; 
clear 

load  hinfcc.mat 

% . End  of  JOSE.M  -  RYC/MGS  % 
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****************************************************************** 


*  This  program  simulates  the  model  with  the  H®  controller  generated* 

•  by  the  modified  Matlab  Robust  Control  Toolbox  program  Jose.m  * 

****************************************************************** 


%  In  order  to  run  this  program  some  of  the  variables  used  in  the 
%  model  simulation  program  are  required.  This  program  simulates  the 
%  model  with  the  H-infinity  controller. 

clg 

format  short 
format  compact 

if  exist(’acp')  =  =  0 
jose 

else 

load  hinfcc.mat 
load  var.mat 
end 

Q  =  diag([f(  1  ),l,f(2),  l,f(3),l,f(4),  l,f(5),  l,f(6),  l,f(7),l,f(8),... 

l,f(9),l,f(10),l].0); 

N  =  200; 


%  Initialize  the  random  inputs  to  the  same  for  each  run. 
rand('normar); 

%  Initialization  of  the  state  and  control  state  vector 

xcl  =  [bg;zeros(bcp)]; 

aida(:,l)  =  xcl(  1:20,1); 

y(l)  =  cg*aida(:,l); 

ww  =  zeros(N,l); 

%  Calculation  of  controller 
ts  =  0.1; 

acl  =  [ag-bg*dcp*cg  -bg*ccp;bcp*cg  acp]; 
bcl  =  [bg:zeros(bcp)]; 

[Phi_cl,Del_cl]  =  c2d(acl,bcl,ts); 


for  n  =  1:N 


%  Plant  simulation  with  the  H-infinitv  controller  of  the  form: 
c/c  xcdot  =  acp’xc  +  bcp*y 

c/c  uc  =  -(ccp*xc  +  dcp*y) 
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xcl(:,n+l)  =  Phi_cl*xcl(:,n)  +  Del_cl*ww(n); 
aida(:,n+  l)=xcl(l:20,n  +  1); 
lccp  =  length(ccp); 

uc  =  -ccp*xcl(2 1 :20  +  lccp,:)-dcp*cg*xcl(  1 :20,:); 

%  Calculation  of  energy’  in  system 
for  mode  =1:10 

E_h(mode,n)  =  0.5*[aida(2*mode-l,n)  aida(2*mode,n)]*... 
Q((2*mode-l):(2*mode),(2*mode-l):(2*mode))*... 
[aida(2*mode-l,n);aida(2*mode,n)]: 
end 
end 

c/c  Plots 
subplot(21 1) 
time  =  0:l:N; 

plot(time,aida(modenbr,:)); 

title([  '  Position  of  Mode  ',num2str(modenbr)]): 

xlabel('seconds');ylabel('Displacement'): 

grid; 

%meta  posl 
pause 

c/c  Plot  of  Total  Dissipation  of  Modes 
for  d  =  1:10 

h_disp(:,d)  =  aida(2*d-l,:)'*bg(2*d); 
end 

:_disp  =  sum(h_disp'); 

%  Plot  of  Total  Displacement 
htime  =  0:l:N;clg: 
plot(htime,t_disp); 

xlabel(Time(sec)  ');ylabel('Amplitude');grid; 

title(Total  Displacement  of  Modes  with  H-infinity  Controller'); 

%meta  hindisp 

pause 

%  Plot  of  Total  Energy 
E_htot  =  sum(E_h); 
time  =  0:ts:19.9; 

subplot(22 1  ),;plot(time,E_htot);grid; 
titIe(Total  Energy'); 

xlabel(Time  (sec)');ylabel('Energy  (in-lbs)'); 
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pause  ;axis; 

9£meta  test2 

%  Plot  of  Total  Energy  per  Mode 
modes*  1:1:10; 

Emh_tot  =  sum(E_h'); 

subplot(222),axis([l  10  0  0.015]);plot(modes,Emh_tot);grid; 
title(’Total  Energy  Per  Mode  ’); 
xlabel('Modes  ');ylabel(’Energy  (in-lbs)');pause 
%meta  energyS 

subplot(223),a\is;plot(uc);grid; 

title(['Control  Input  for  Mode  ',num2str(modenbr)]); 

xlabel('  Time  (sec)');ylabel('Energy  (in-lbs)'); 

%meta  control  1 
pause 

subplot(224), 

plot(diag(uc'*uc));title(’Control  Energy');grid; 
xlabel(Time  (sec)');ylabel('Energy  (in-lbs)'); 

%meta  ful!w2 
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******************************************************** 


*  LQGCTCR.M 

*  This  program  is  designed  to  generate  a  reduced  order  LQ  regulator  * 

*  controller.  The  number  of  modes  to  be  controlled  is  determined  when  * 

*  the  user  is  queried  by  the  program.  * 

************************************************************************ 

%  In  order  to  run  this  program  some  of  the  variables  used  in  the 
%  model  simulation  program  are  required. 

if  exist(’modenbr')  =  =  1 
load  var.mat 
else 
var 
end 


%Plant  matrices  for  20th  order  system  for  full  order 
%state  feedback  controller 

dv  =  zeros(  19,1); 
for  i  =  1:19 
if(rem(i,2)'  =0); 

dv(i,l)=  1; 
end 
end 


A  =  diag(dv,  1 )  +  diag([-0. 1048,0,-0. 1206,0,-0. 1435,0,-0. 1452,0,... 

-0.15697, 0,-0.1643, 0,-0.1662, 0,-0.1724, 0,-0.1864,0, -0.1972], -1)  +  ... 
diag([0,-0.00324, 0,-0.00347,0,-0.00378, 0,-0.00381,0,-0.00396,0,... 

-0.00405, 0,-0.00408, 0,-0.00415, 0,-0.00432,0, -0.00444], 0); 

B  =  [0;b(  l);0;b(2);0;b(3);0;b(4);0;b(5);0;b(6);0;b(7);0;b(8);0;b(9);... 

0;b(  10)]; 

Q  =  diag([f(l),l,f(2),l,f(3),l,f(4),l,f(5),l,f(6),l,f(7),l,f(8),... 
l,f(9),l,f(10),l],0); 

R=  l; 

%  Number  of  modes  to  be  controlled 

nbr_mode  =  input(’How  many  modes  do  you  desire  to  control  (1:10)?  ’); 

%  The  following  procedure  is  for  a  reduced  order  controller. 

%  The  user  selects  the  number  of  modes  to  be  controlled. 

ts  —  1;  %Sampling  Time 
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[Phi,Del]  =  c2d(A,B,ts); 

aida(:,l)  =  B;  %  Initialization 

Phi_r  =  Phi(l:nbr_mode*2,l:nbr_mode*2); 

Del_r  =  Del(l:nbr_mode*2); 

0_r  =  0(1 :  nbr_mod  e  *2,1 :  nbr_mod  e  *  2 ) ; 

K_r  =  dlqr(Phi_r,Del_r,0_r,R); 

K  =  [K_r  zeros(l,20-nbr_mode*2)]; 

for  n=  1:100  %  Time 

u  =  -K*aida(:,n);  c/c  No  control  added 

for  mode  =  1:10 

Phimode  =  Phi(2*mode-l:2*mode,2*mode-l:2*mode); 

Deimode  =  Del(2*mode-l:2*mode); 

Oe  =  Q(2*mode-l:2*mode,2*mode-l:2*mode);  %  Energy  O  matrix 

%  Aida  gives  the  matrix  for  twenty  states  for  each  point  in  time. 

%  Next  actual  state  calculated 

aida(2*mode-l:2*mode,n+  l)  =  Phimode*aida(2*mode-l:2*mode,n)... 
f  Deimode  *  u  ;  %  control  added 

c/c  At  each  point  in  time  the  energy  of  each  mode  is  calculated  using 
%  Energy  in  specified  mode  =  0.5[aida  aida_d]Q[aida 
c/c  aida_d] 

E(mode,n)  =  0.5*[aida(2*mode-l,n)  aida(2*mode,n)]  *... 

Oe  *[aida(2*mode-l,n);aida(2*mode.n)]; 


end 

end 

%  Plot  of  Total  Dissipation  of  Modes 
for  m  =  1:10 

out_d(:,m)  =  aida(2*m-l,:)'*B(2*m); 
end 

outtot  =  sum(outd’); 

%  Plot  of  Total  Displacement 
time_o  =  0:l:100; 
plot(time_o,out_tot); 

xlabel(Time(sec)  ');ylabel('Amplitude');grid; 
title(Total  Displacement  of  Modes’); 
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gtext(['R  =’,  num2str(R)  ]); 
meta  redd 
clg; 


E_total  =  sum(E); 
time  =  0:1:99; 

subplot(211),plot(time,E_total);grid; 
title(Total  Energy  with  Control'); 
xlabel(Time'); 
gtext(['R  =  \num2str(R)]); 

gtext(['Number  of  Modes  Controlled:  ',num2str(nbr_mode)]); 
pause 

%  Plot  of  Total  Energy  per  Mode 
Em_totaI  =  sum(E'); 
modes=  1:1:10; 

subplot(212),plot(modes,Em_total);grid; 
title(Total  Energy  per  Mode  with  Control'); 
gtext(['R  =  ’,num2str(R)]); 
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•ft********************************************************************** 


*  KALMANR.M 

*  This  program  generates  the  reduced  order  controller  for  the  Kalman  * 

*  filter.  * 


************************************************************************ 


9c  In  order  to  run  this  program  some  of  the  variables  used  in  the 
9c  model  simulation  program  are  required. 

9c  This  program  simulates  the  model  with  estimated  states  using 
9c  steady  state  Kalman  filter. 

clg 

format  short  e 
format  compact 


if  exist('b')  =  =  1 
load  var.mat 
else 
var 
end 

9c  Plant  matrices 
dv  =  zeros(19,l); 
for  i  =  1:19 
if(rem(i,2)‘  =0); 

dv(i,l)=  1; 
end 
end 

A  =  diag(dv,l)  +  diag([-0.1048,0,-0.1206, 0,-0.1435,0, -0.1452,0,... 

-0. 1 5697,0,-0. 1643,0,-0. 1662,0,-0. 1724,0.-0  1 864,0,-0. 1 972].- 1  )  +  ... 
diag([0,-0.003 24,0,-0.00347,0, -0.00378,0,-0.0038 1,0,-0.00396,0.... 

-0.00405,0,-0.00408,0,-0.0041 5, 0,-0.00432,0,-0.00444],0); 

B  =  [0;b(l);0;b(2);0;b(3);0;b(4);0,b(5);0;bf6);0;b(7);0;b(8);0;b(9);... 
0;b(  10)]; 

C  =  [b(l)  0  b(2)  0  b(3)  0  b(4)  0  b(5)  0  bVo)  0  b(7)  0  ... 

b(8)  0  b(9)  0  b(  10)  0]; 

D  =  0; 

Q  =  diag([f(l),l,f(2).l,f(3),l,f(4),l,f(5),l,f(6),l,f(7),l,f(8),... 

1  ,f(9).  1  ,f(  1 0).  1  ],0); 

R=  1; 

N  =  200; 

nbr_mode  =  input(’Number  of  modes  to  control:  '); 
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%  Initialize  the  random  inputs  to  the  same  for  each  run. 
rand('normar); 

rand(’seed',0);  9c  Sets  the  seed  to  0  when  Matlab  is  entered 

%  Initialization  of  the  state  and  the  estimate,  x_hat 

%  and  the  error  covariance  matrix. 

aida  hat(:,l)  =  zeros(B); 

aida(:,l)  =  zeros(B); 

old_P  =  zeros(A); 

9c  Desired  state 
aida_d  =  zeros(B); 

ts  =  l;  %Sampling  Time 

W  =  40; 

V  =  0.005; 

[Phi, Del]  =  c2d(A,B,ts); 

G  =  real(dlqe(Phi,Del,C,W,V));  %  Steady  State  Optimal  Gain 

K  =  real(dlqr(Phi,Del,Q,R));  9c  Gain  matrix  K 

%For  reduced  order  matrices 

Phi_kr  =  Phi(l:nbr_mode*2,l:nbr_mode*2); 

Del_kr  =  Del(l:nbr_mode*2); 

Q_kr  =  Q(l:nbr  mode*2,l:nbr_mode*2); 

K_kr  =  dlqr(Phi_kr,Del_kr,Q_kr,R); 

K  =  [K_kr  zeros( l,20-nbr_mode*2)]; 


for  k=  1:N 


%  Plant  simulation 

ww(k)  =  zeros(N,  l);ww(l)  =  1; 

u(k)  =  -K*(aida_hat(:,k)-aida_d); 

aida(:,k+ l)  =  Phi*aida(:,k)  +  Del*u(k)  +  Del*ww(k); 

w(k+  l)  =  sqrt(V)*rand; 

y(k+ l)  =  C*aida(:,k+ 1)  +  w(k+l); 

%  Estimates  updated 

aida_hat(:,k+ 1)  =  Phi*aida_hat(:,k)  +  Del  *u(k); 
y_hat(k+l)  =  C  *aida_hat(:,k+ 1); 

aida_hat(:,k+ 1)  =  aida_hat(:,k+ 1)  +  G*(y(k+ l)-y_hat(k+ 1)); 
for  mode  =  1:10 

E_k(mode,k)  =  0.5 *[aida(2* mode- l,k)  aida(2*mode,k)],t... 
0(2*  mode- 1:2*  mode,  2*  mode- 1 :2*  mode)*... 
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[aida(2*mode-l,k);aida(2*mode,k)]; 


end 

end 


u(N+  l)  =  -K*aida_hat(:,N  + 1); 
t(N+  1)  =  N+  1; 

%  Plots 
subplot(211) 

plot(t,aida(modenbr,:).’o’,t,aida_hat(modenbr, :),'*'); 
title(['Estimated  and  Actual  Position  of  Mode  ',num2str(modenbr)]); 
xlabel(Time(sec)');ylabel('Displacement'); 
grid;f/fgtext(['W  =  ',num2str(W)]);gtext(['V  =  ’,num2str(V)]); 

%gtext('w/  rand  plant  disturbance  and  meas  noise’); 

%gtext([’Reduced  Order  Controller/Kalman  Filter;', num2str(nbr_mode)]); 
pause 

c/c  Plot  of  Total  Dissipation  of  Modes 
for  d  = 1:10 

k_disp(:,d)  =  aida(2*d-l,:)'*B(2*d); 
kh_disp(:,d)  =  aida_hat(2*d-l,:)'*B(2*d); 
end 

tdisp  =  sum(k_disp’); 
thdisp  =  sum(kh_disp'); 

%  Plot  of  Total  Displacement 
ktime  =  0:l:N;subplot(212) 
plot(ktime,t_disp,'--',ktime,th_disp, '-.'); 
xlabel(Time(sec)  ');ylabel('Amplitude’);grid; 

title(’  Total  Displacement  of  Modes  with  Reduced  Order  Kalman  Filter'); 

meta  kalrdiw3 

pause 

clg 

K_energv  =  0.5*(aida’*Q*aida); 

Khenergy  =  0.5*(aida_hat'*Q'aida_hat);subplot(21 1) 

plot(ktime,diag(K_energy),'o',ktime,diag(Khenergy),’*');grid; 

title('Actual(o)  and  Est.(*)  Energy  w/  Reduced  Order  Kalman  Filter’); 

xlabel(Time’);ylabe!(’Magnitude'); 

pause; 


%  Plot  of  Total  Energy  per  Mode 
modes  =  1:1:10; 
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Emk_tot  =  sum(E_k'); 
subplot(2l2),plot(modes,Emk_tot);grid; 

title(’Total  Energy  Per  Mode  with  Reduced  Order  Kalman  Filter'); 
xlabel('Modes  ');ylabel(’Amplitude');pause 
meta  kalrepw3 
clg 

plot(u);grid;xlabel('Time  (sec)');ylabel('Magnitude'); 
title('Control  Input  with  Reduced  Order  Kalman  Filter');pause 
clg 

plot(t,diag(u'*u)); 

grid;title('Control  Energy  with  Reduced  Order  Kalman  Filter'); 
xlabel(Time  (sec)');ylabel('Amplitude'); 
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************************************************************************ 


»  JOSER.M 

*  This  program  generates  the  reduced  order  H°°  controller  for  * 

*  the  Hinfcr.m  simulation  program.  * 

************************************************.************************ 


%  JOSER  Large  space  structure  design  demonstration 
c/c 

%  R.  Y.  Chiang  &  M.  G.  Safonov  3/88 
%  Copyright  (c)  1988  by  the  Math  Works,  Inc. 

%  All  Rights  Reserved. 


% . 

clc 

disp('  <  <  Demo  #3:  MIMO  Large  Space  Structure  Design  Example  >  >') 

dispC  _ ') 

disp('  Secondary  Mirror  - — >  -ooo-  *') 

dispC  f\/\  I') 

dispf  I  X  |  |’) 

dispC  I  /  \  I  I') 

dispC  .  I') 

dispC  I  \  /  I  I') 

dispC  |  X  |  !') 

dispC  I  /  \  I  ’) 

dispC  Lens . >  — O—  7.4  Meters’) 

dispC  |  |  ’) 

dispC  I  \  /  !  D 

dispC  I  \  /  I  I') 

dispC  I  \  /  t  C) 

dispC  I  \  /  I  I') 

dispC  |  X  |  0 

dispC  I  /  \  I  I') 

dispC  I  /  \  I  I') 

dispC  /  /  \  \  I') 

dispC  Primary  Mirror  -->  (_000000 _ \ _ )  |') 

dispC  \  /I') 

dispC  \ _ / _ v _ ') 

dispC  ’) 

dispC  (strike  a  key  to  continue  ...)’) 

pause 

if  exist(  b  )  =  =  1 


load  var.mat 
else 
var 
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end 


nbr_mode  =  input('Enter  the  number  of  modes  to  control:  '); 

format  short 
dv  =  zeros(19,l); 
for  i  =  1:19 

if  (rem(i,2)*=0); 

dv(i,l)=  1; 
end 
end 


ag  =  diag(dv,l)  +  diag([-0.1048, 0,-0.1206, 0,-0. 1435,0,-0. 1452,0,... 

-0. 15697,0,-0. 1 643,0,-0. 1662,0,-0. 1724,0,-0. 1 864,0,-0. 1972],- 1 )  + ... 
diag([0, -0.00324, 0,-0.00347, 0,-0.00378,0, -0.0038 1,0, -0.00396,0,... 

-0.00405, 0,-0.00408, 0,-0.00415,0, -0.00432, 0,-0.00444],0): 
bg  =  [0;b(  l);0;b(2);0;b(3);0;b(4);0;b(5);0;b(6);0;b(7);0;b(8);0;... 
b(9);0;b(  10)]; 

eg  =  [b(l),0,b(2),0,b(3),0,b(4),0,b(5),0,b(6),0,b(7),0,b(8),0.... 

b(9),0,b(10),0]; 
dg  =  0; 
clc 

%  To  reduce  the  order  of  the  model,  the  size  of  the  matrix  is 

%  determined  the  number  of  modes  to  control 

ag_r  =  ag(  1  :nbr_mode  *2, 1  :nbr_mode  *2); 

bg_r  =  bg(l:nbr_mode*2); 

cg_r  =  cg(l,l:nbr_mode*2); 

dg_r  =  0; 

dispC  ') 

dispC  State-Space  of  the  Large  Space  Structure:') 

disp('  (after  colocated  rate  feedback  and  model  reduction)’) 

%ag 

%bg 

%disp('  (strike  a  key  to  continue  ...)') 

%pause 

%clc 

%cg 

%dg 

%disp('  (strike  a  key  to  continue  ...)') 

%pause 

%clc 

dispC  ’) 

dispC  ’) 

dispC  Poles  of  the  Plant  (stable):') 
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dispC  ') 

dispC  . - . . . ’) 

dispC  poleg  =  eig(ag)  c/c  Computing  the  poles  of  the  plant') 

dispC  . - . — *) 

poleg  =  eig(ag_r) 
dispC  ') 
dispC  ') 

dispC  (strike  a  key  to  continue  ...)') 

pause 

clc 

dispC  ') 


dispC  ') 

dispC  Transmission  Zeros  of  the  Plant  (minimum  phase):') 
dispC  ’) 


dispC  . . ’) 

dispC  tzerog  =  tzero(ag,bg,cg.dg)  c/c  Computing  the  transmission  zeros') 
dispC  . ') 

tzerog  =  tzero(ag_r,bg_r,cg_r,dg)  c/c  Computing  the  transmission  zeros 


dispC  ') 
dispC  ’) 
dispC  ') 

dispC  (strike  a  key  to  continue  ...)') 

pause 

clc 

dispC  ') 

dispC  -  -  -  Computing  SV  Bode  plot  of  the  open  loop  plant  -  -  -') 
w  =  logspace(-3,5,100); 

svg  =  sigma(ag_r,bg_r,cg_r,dg.l,w);  svg  =  20*logl0(svg): 

dispC  ') 

dispC  ') 

dispC  ') 

dispC  ') 

dispC  ') 

dispC  (strike  a  key  to  see  the  SV  Bode  plot  of  G(s)  ...)') 

pause 

semilogx(w.svg) 

title('MIMO  Larse  Space  Structure  Open  Loop') 
xlabel('Frequency  -  Rad/Sec') 
ylabe!(’SV  -  db’)' 
grid 

pause;%meta  jackl 
clc 

dispC  ') 

dispC  <  <  Design  Specifications  >  >  ') 
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dispC  ') 

dispC  1).  Robustness  Spec.  :  closed  loop  bandwidth  --  200  r/s  (30  hz)’) 
disp(’  Associated  Weighting:') 

disp('  ') 

dispC  -1  200') 

dispC  W3(s)  =  . *  I  (fixd)') 

dispC  s  2x2') 

dispC  ') 
dispC  ') 

dispC  2).  Performance  Spec.:  sensitivity  reduction  of  at  least  10:1') 
disp('  UP  to  approx.  10  r/s') 

dispC  Associated  Weighting:') 

dispC  ') 

dispC  -1  -1  0.01  (1  +  s/lOrr) 

dispC  Wl(s)  =  Gam  *  .  *  I') 

dispC  (1  +  s/500)  *2  2x2’) 

dispC  ') 

dispC  where  "Gam"  in  our  design  goes  from  1  -->  1.5.') 

coef=input('Enter  the  coef  value'); 
fac  =  500/coef; 
k  =  200/fac; 

nuw3i  =  [0  k];  dnw3i  =  [1  0]; 

nuwli  =  conv([l/(  10/fac)  1],[  l/(  10/fac)  1]); 

dnwli  =  100*conv([l/coef  l],[l/coef  1]); 

%k  =  200; 

%nuw3i  =  [0  k];dnw3i  =  (1  0]; 

svw3i  =  bode(nuw3i,dnw3i,w);  svw3i  =  20*logl0(svw3i); 

%nuwli  =  conv([l/10  1].[1/10  l]);dnwli  =  1 00 *conv([  1/500  1],[  1/500  1]); 
svwli  -  bode(nuwli,dnwli,w);  svwli  =  20*logl0(svwli); 
dispC  ') 
dispC  ') 

dispC  (strike  a  key  to  see  the  plot  of  the  weightings  ...)') 

pause 

axis([0  5  -40  40]) 

semilogx(w,svwli,w,svw3i) 

grid 

title('MIMO  LSS  Design  Example  -  Design  Specifications') 

xlabel('Frequency  -  Rad/Sec') 

ylabel('l/Wl  &  1/W3  -  db') 

text(2,-20,'Sensitivity  Spec.-  1/Wl(s)') 

text(2, 10, 'Robustness  Spec.-  1/W3(s)') 

pause;%meta  jack  1 

axis 

clc 
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dispC 
dispO  ') 

dispC 

dispC  ') 

dispC 

dispC  ’) 

dispC 

dispC  ’) 

dispC 

dispC 

dispC  ') 

dispC 

dispC 

dispC  ') 

dispC 

dispC 

dispC 

dispC 

dispC 

dispC  ') 

dispC  ') 

dispC  ') 

dispC 

pause 

clc 

dispC  ’) 

dispC  ') 

dispC 

dispC  ') 

dispC 

dispC 

dispC 

dispC 

dispC 

dispC 

dispC 

dispC 

dispC 

dispC 
dispC 
dispC 
dispC  ’) 
dispC  ') 
dispC 


<  <  Problem  Formulation  >  >') 


Form  an  augmented  plant  P(s)  with  these  two  weighting  function*:') 

1) .  W1  penalizing  error  signal  "e"') 

2) .  W3  penalizing  plant  output  V") 


and  find  a  stabilizing  controller  F(s)  such  that  the  Hinf-norm') 
of  TF  Tylul  is  minimized  and  less  than  one,  i.e.') 


min  |Tylul  |  <  1.’) 

F(s)  inf) 


where  ') 

Tv  lu 


-i|') 

|  Gam*\VT*(I  +  GF)  |  =  |  Gam  *  \V1  *  S  |') 
-1|  I  W3  *  (I  -  S)  D 
\V3*GF*(I  +  GF)  !') 


(strike  a  key  to  continue  ...)’) 


<  <  DRSIGN  PROCEDURE  >  >') 

************************,*,  »*,^ 

[Step  1].  Do  plant  augmentation  (run  augtf.m  or  *’) 
augss.m)  *’) 

•') 

[Step  2].  Balanced  the  augmented  plant  for  better  *') 

numerical  condition  (run  OBALREAL.M)  *') 

*') 

[Step  3].  Do  H-inf  synthesis  (run  HINF.M)  *') 

*') 

[Step  4j.  Redo  the  plant  augmentation  and  balancing  *') 

for  a  new  "Gam"  -->  1.5  and  rerun  HINF.M  *') 

**************************** 


(strike  a  key  to  continue  ...)') 
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pause 

clc 

dispC  ') 

dispC  ') 

dispC  Assign  the  cost  coefficient  "Gam"  -->  1  ’) 

dispC  ’) 

dispC  this  will  serve  as  the  baseline  design 

dispC  •) 
dispC  ’) 

Gam  =  input('  Input  the  cost  coefficient  "Gam"  =  ’); 


dispC  ') 

dispC  . ') 

dispC  c/c  Plant  augmentation  of  the  LSS:') 

dispC  sysg  =  [ag  bg;cg  dgl;  xg  =  4;’) 

dispC  wl  =  [Gam*dnwli;nuwli;Gam*dnwii;nuwli];') 

dispC  w2  =  []:) 

dispC  w3  =  [dnw3i;nuw3i;dnw3i;nuw3i]:') 

dispC  [A,B1,B2,C1,C2,D11,D12,D21,D22]  =  augtf(sysg,xg,wl,w2,w3);') 

dispC  . ') 

sysg  =  [ag_r  bg_r;cg_r  dg];  xg  =  nbr_mode*2; 

wl  =  [Gam*dnwli;nuwli];  %Gam*dnwli;nuwli]; 

w2  =  []; 


w3  =  [dnw3i;nuw3i];  %  dnw3i;nuw3i]: 

[A,B1,B2,C1,C2,D1 1,D12,D21,D22]  =augtf(svsg,xg,wl,w2,w3); 
dispC  ’) 


dispC  -  -  -  State-Space  (A.B1,B2,C1,C2,D11,D12,D21,D22)  is  ready  for  ) 
dispC  the  Small-Gain  problem  -  -  -') 

dispC  ’) 

dispC  . ') 

dispC  [aa,bb,cc,mm,tt]  =  obalreal(A[Bl  B2],[C1;C2])  c/c  Balancing  P(s)') 
dispC  A  =  aa;  B1  =  bb(:,l:2);  B2  =  bb(:,3:4);  ') 
dispC  Cl  =  cc(  1:4,:);  C2  =  cc(5:6,:):') 
dispC  . --') 


%[aa,bb,cc,mm,tt]  =  obalreal(A,[Bl  B2].[Cl;C2j); 

%A  =  aa;  B1  =  bb(:,l:2);  B2  =  bb(:,3:4);  Cl  =  cc(l:4,:);  C2  =  cc(5;6,;); 
dispC  ') 
dispC  ') 

dispC  (strike  a  key  to  continue  ...)') 

pause 

clc 

dispC  ’) 
dispC  ’) 
dispf  ’) 

dispC  . - . ') 
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disp('  hinf  %  Running  script  file  HINF.M  for  H-inf  optimization’) 

dispC  . - . ') 

hinfr 
dispC  ') 
dispC  ') 

disp('  (strike  a  key  to  continue  ...)') 

pause 

pltoptr  c/c  Preparing  singular  values  for  plotting 

svwlil  =  svwli;  hsvsl  =  svs;  hsvtl  =  svt;  hsvttl  =  svtt; 
dispC  ') 
dispC  ') 

dispC  (strike  a  key  to  continue  ...)') 


pause 


clc 

dispC  ') 
dispC  ’) 

dispC  After  a  few  iterations,  we  found  a  new  Gam  of  1.5  can  push  the’) 
dispC  ') 

dispC  H-inf  cost  function  close  to  its  limit.  ') 
dispC  ') 
dispC  ') 

dispC  Input  "Gam"  -->  1.5,  and  try  HINF  again . ’) 

dispC  ') 

dispC  ') 

Gam  =  input('  Input  the  cost  coefficient  "Gam"  =  '); 

dispC  ') 

dispC  . ’) 

dispC  %  Adjust  plant  augmentation:’) 

dispC  wl  =  [Gam*dnwli;nuwli;Gam*dnwli;nuwli];') 

disp('  [A,B1,B2,C1,C2,D1 1,D12,D21,D22]  =  augtf(sysg,xg.wl,w2,w3);') 

dispC  . - . ') 

wl  =  [Gam*dnwli;nuwli];  %Gam*dnwli;nuwli]; 


[A,B1.B2,C1,C2.D1 1,D12,D21.D22]  =  augtf(sysg,xg.wl,w2,w3); 
%[aa,bb,cc,mm,tt]  =  obalreal(A,[Bl  B2],[C1;C2]); 

%A  =  aa;  B1  =  bb(:,l:2);  B2  =  bb(:,3:4);  Cl  =  cc(l:4,:);  C2  =  cc(5:6,:); 
dispC  ') 
dispC  ') 

dispC  (strike  a  key  to  continue  ...)') 

pause 

hinf 

dispC  ') 
dispC  ') 

dispC  (strike  a  key  to  continue  ...)') 

pause 
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pltoptr 

svwli2  =  svwli;  hsvs2  =  svs;  hsvt2  =  svt;  hsvtt2  =  svtt; 
dispC  ') 
disp('  ') 

disp('  (strike  a  key  to  see  the  plots  of  the  comparison  ...)') 

pause 

semilogx(w,svwlil,w,hsvsl,w,svwli2,w,hsvs2) 
title('H-inf  LSS  Design  --  1/W1  &  Sensitivity  Func.') 
xlabel(’Frequencv  -  Rad/Sec') 
ylabel('SV  -  db’)' 
grid; 

text(0.002,0, 'H-inf  (Gam  =  1)  — >  H-inf  (Gam  =  1.5)') 

pause;%meta  jackl 

se  mil  ogx(  w,svw3  i,w,  hsvt  1  ,w,  hsvt2 ) 

title(’H-inf  LSS  Design  --  1/W3  &  Comp.  Sens.  Func.’) 

xlabel(’Frequencv  -  Rad/Sec') 

ylabelfSV  -  db')' 

grid 

text(0.002,0,'H-inf  (Gam  =  1)  — >  H-inf  (Gam  =  1.5)') 

pause;%meta  jackl 

semilogx(w,hsvttl,w,hsvtt2) 

title('H-inf  LSS  Design  -  Cost  function  Tylul') 

xlabel('Frequency  -  Rad/Sec') 

ylabel('SV  -  db’) 

grid 

text(0.002, -10, 'H-inf  (Gam  =  1)  — >  H-inf  (Gam  =  1.5)') 
pause;%meta  jackl 

clc 

disp('  ’) 
disp('  ') 

disp(’  <<  8-State  H-inf  Controller  (Gam  =  1.5)  >>') 

dispC  ') 

disp('  Poles  of  Controller  :') 
polecp  =  eig(acpr) 
disp(’  ') 

disp(’  (strike  a  key  to  continue  ...)') 

pause 

clc 

disp('  State-Space  of  the  8-State  H-inf  Controller:’) 
disp('  First  6  columns  of  the  A  matrix:') 
acp_r(:,l:nbr_mode*2) 

dispC  (stiike  a  key  to  continue  ...)’) 

pause 
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clc 

dispf  ') 

%disp('  Last  two  columns  of  the  A  matrix:') 

%acp(:,7:8) 

disp(’  (strike  a  key  to  continue  ...)') 

pause 

clc 

bcp_r 

dispC  (strike  a  key  to  continue  ...)') 

pause 

clc 

ccp_r 

dcp_r 

dispC  (strike  a  key  to  continue  ...)') 

pause 

clc 

dispC  ’) 

dispC  Poles  of  closed-loop  TF  matrix  Tylul:') 
poletyu  =  eig(acl) 

dispC  (strike  a  key  to  continue  ...)') 

pause 

Cc 

% . End  of  JOSEDEMO.M  -  RYC/MGS  % 

save  hinfcr.mat  acp_r  bcp_r  ccp_r  dcp  r  b  ag  bg  eg  dg  nbr_mode  modenbr  f; 
clear 

load  hinfcr.mat 
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************************************************************************** 


*  HINFCR.M 

*  This  program  simulates  the  model  with  the  reduced  order  H*>  controller  * 

%  In  order  to  run  this  program  some  of  the  variables  used  in  the 
%  model  simulation  program  are  required.  This  program  simulates  the 
%  model  with  the  Reduced  Order  H-infinity  controller. 

clg 

format  short  e 
format  compact 

if  exist('acp’)  =  =  0 
joser 
else 

load  hinfcr.mat 
load  var.mat 
end 

%  Plant  matrices 
dv  =  zeros(19,l); 
for  i  =  1:19 
if(rem(i,2)*=0); 

dv(i,l)=  1; 
end 
end 

Q  =  diag([f(l),l,f(2),l,f(3),l,f(4),l,f(5),l,f(6).l,f(7),l,f(8),... 

l,f(9),  l,f(  10),  1],0); 

N  =  200; 


%  Initialize  the  random  inputs  to  the  same  for  each  run. 
rand('normal'); 

rand('seed',0);  %  Sets  the  seed  to  0  when  Matlab  is  entered 

%  Initialization  of  the  state  and  control  state  vector 

aida(:,l)  =  bg; 

xcl_r  =  zeros(  [bg;bcp] ) ; 

y(l)  =  cg*aida(:,l); 

ww  =  zeros(N,l);ww(l)  =  1; 

%Calculation  of  controller  and  Conversion  to  discrete  time 
ts  =  0. 1 ; 

acl_r  =  [ag-bg*dcp*cg  -bg*ccp;bcp*cg  acp]; 
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bcI_r  =  [bg;zeros(bcp)]; 

[Phi_clr,Del_clr]  =  c2d(acl_r,bc!_r.ts); 

for  n=  1:N 

c/c  Plant  simulation  with  the  H-infinity  controller  of  the  form: 

9  c  xc_dot  =  acp*xc  +  bcp*y 

9c  uc  =  ccp*xc  +  dcp*y 

%  Controller  state  updated  and  input  updated 
xc!_r(:,u+l)  =  Phi_clr*xcl_r(:.n)  +  Del_clr*ww(n); 
aida(:,n  +  l)  =  xcl_r(l:20,n  +  1); 
lccp_r  =  length(ccp); 

uc_r  =  -ccp*xcl_r(2 1 :20  +  lccp_r,:)-dcp*cg*xcl_r(  1 :20,:); 

9 c  Calculation  of  energy-  in  system 
for  mode  =  1:10 

T_hr(mode.n)  =  0.5  *[aida(2*  mode- l.n)  aida(2*mode,n)]*... 
0(2*  mode- 1:2  ‘mode.  2*  mode- 1:2*  mode)*... 

[aida(2*  mode- l.n);aida(2*  modern)]; 
end 

end 

9c  Plots 
time  =  0:1:N: 

plot(time.aida(modenbr,:)); 

title( [  Position  of  Mode  ',num2str(modenbr)]): 

xlabel('seconds');yIabeI('DispIacement'); 

grid;9cgtext('w/  impulse  plant  disturbance  ’);9cmeta  jack 

pause 

9c  Plot  of  Total  Dissipation  of  Modes 
for  d  =  1:10 

h_disp(:,d)  =  aida(2*d-l,:)'*bg(2*d); 
end 

tdisp  =  sum(h_disj:  ’); 

9c  Plot  of  Total  Displacement 
htime  =  0:  t:N;clg; 
plot(htime,t_di: ,  ); 

xlabelCTime  (sec)  ');y!abel('Amplitude');grid: 

title(Total  Displacement  of  Modes  w/  H-infinity  Reduced  Controller' 
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%meta  hinfcr 
pause 

%  Plot  of  Total  Energy 
time  =  0:ts:19.9; 

subplot(221),plot(time,sum(E_hr)); 
title(Total  Energy  ');grid; 
xlabel(Time  (sec)');ylabel('Energy  (in-lbs)'); 
pause; 

%meta  testr2 

%  Plot  of  Total  Energy  per  Mode 
modes  =  1:1:10; 

Emh_tot  =  sum(E_hr’);subplot(222) 

subplot(222),axis([l  10  0  0.015]);plot(modes,Emh_tot);grid; 
title(Total  Energy  Per  Mode  '); 
xlabel('Modes  ');ylabel('Energy  (in-lbs)');pause 
%meta  hinfcr 

subplot(224), 

plot(diag(uc_r'*uc_r));title('Control  Energy');grid; 
xlabel('  Time  (sec)');ylabel('  Energy  (in-lbs)'); 
pause 

meta  redl5_7; 

clg;plot(uc_r);grid; 
title('Control  Input'); 

xlabelO  Time  (sec)');ylabel('Energy  (in-lbs)'); 

%meta  contred7 
pause 
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